107 auto outputBlock = std::make_unique<QgsRasterBlock>();
111 return outputBlock.release();
114 std::shared_ptr< QgsRasterBlock > inputBlock(
mInput->block( mBand,
extent, width, height, feedback ) );
116 if ( !inputBlock || inputBlock->isEmpty() )
119 return outputBlock.release();
122 std::shared_ptr< QgsRasterBlock > alphaBlock;
127 if ( !alphaBlock || alphaBlock->isEmpty() )
130 return outputBlock.release();
135 alphaBlock = inputBlock;
140 return outputBlock.release();
143 if ( width == 0 || height == 0 )
144 return outputBlock.release();
149 float cellXSize =
static_cast<float>(
extent.width() ) / width;
150 float cellYSize =
static_cast<float>(
extent.height() ) / height;
151 float zenithRad =
static_cast<float>( std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0 );
152 float azimuthRad =
static_cast<float>( -1 * mLightAzimuth * M_PI / 180.0 );
153 float cosZenithRad = std::cos( zenithRad );
154 float sinZenithRad = std::sin( zenithRad );
157 float cos_alt_mul_z = cosZenithRad *
static_cast<float>( mZFactor );
158 float cos_az_mul_cos_alt_mul_z = std::cos( azimuthRad ) * cos_alt_mul_z;
159 float sin_az_mul_cos_alt_mul_z = std::sin( azimuthRad ) * cos_alt_mul_z;
160 float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0f * cos_az_mul_cos_alt_mul_z;
161 float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0f * sin_az_mul_cos_alt_mul_z;
162 float square_z =
static_cast<float>( mZFactor * mZFactor );
163 float sin_altRadians_mul_254 = 254.0f * sinZenithRad;
166 float sin_altRadians_mul_127 = 127.0f * sinZenithRad;
168 float cos225_az_mul_cos_alt_mul_z_mul_127 = -32.87001872802012f * cos_alt_mul_z;
169 float cos_alt_mul_z_mul_127 = 127.0f * cos_alt_mul_z;
182 if ( source.isEmpty() )
190 std::chrono::time_point<std::chrono::system_clock> startTime( std::chrono::system_clock::now() );
197 std::size_t inputDataTypeSize = inputBlock->dataTypeSize();
198 std::size_t outputDataTypeSize = outputBlock->dataTypeSize();
201 switch ( inputBlock->dataType() )
204 typeName = u
"unsigned char"_s;
207 typeName = u
"unsigned int"_s;
210 typeName = u
"short"_s;
213 typeName = u
"unsigned int"_s;
219 typeName = u
"float"_s;
222 throw QgsException( u
"Unsupported data type for OpenCL processing."_s );
227 source.replace(
"__global float *scanLine"_L1, u
"__global %1 *scanLine"_s.arg( typeName ) );
231 std::size_t scanLineWidth( inputBlock->width() + 2 );
232 std::size_t inputSize( inputDataTypeSize * inputBlock->width() );
235 std::size_t bufferWidth( width + 2 );
236 std::size_t bufferSize( inputDataTypeSize * bufferWidth );
241 auto scanLine = std::make_unique<QgsRasterBlock>( inputBlock->dataType(), scanLineWidth, 1 );
248 std::vector<float> rasterParams;
249 rasterParams.push_back( inputBlock->noDataValue() );
250 rasterParams.push_back( outputBlock->noDataValue() );
251 rasterParams.push_back( mZFactor );
252 rasterParams.push_back( cellXSize );
253 rasterParams.push_back( cellYSize );
254 rasterParams.push_back(
static_cast<float>(
mOpacity ) );
257 rasterParams.push_back( cos_az_mul_cos_alt_mul_z_mul_254 );
258 rasterParams.push_back( sin_az_mul_cos_alt_mul_z_mul_254 );
259 rasterParams.push_back( square_z );
260 rasterParams.push_back( sin_altRadians_mul_254 );
263 rasterParams.push_back( sin_altRadians_mul_127 );
264 rasterParams.push_back( cos225_az_mul_cos_alt_mul_z_mul_127 );
265 rasterParams.push_back( cos_alt_mul_z_mul_127 );
268 rasterParams.push_back(
static_cast<float>( qBlue( defaultNodataColor ) ) );
269 rasterParams.push_back(
static_cast<float>( qGreen( defaultNodataColor ) ) );
270 rasterParams.push_back(
static_cast<float>( qRed( defaultNodataColor ) ) );
271 rasterParams.push_back(
static_cast<float>( qAlpha( defaultNodataColor ) ) / 255.0f );
274 rasterParams.push_back(
static_cast<float>( mMultiDirectional ) );
276 cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(),
true,
false,
nullptr );
277 cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize,
nullptr,
nullptr );
278 cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize,
nullptr,
nullptr );
279 cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize,
nullptr,
nullptr );
280 cl::Buffer *scanLineBuffer[3] = { &scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer };
282 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, outputDataTypeSize * width,
nullptr,
nullptr );
284 static std::map<Qgis::DataType, cl::Program> programCache;
285 cl::Program program = programCache[inputBlock->dataType()];
286 if ( !program.get() )
290 program = programCache[inputBlock->dataType()];
297 auto kernel = cl::KernelFunctor< cl::Buffer &, cl::Buffer &, cl::Buffer &, cl::Buffer &, cl::Buffer & >( program,
"processNineCellWindow" );
301 std::vector<int> rowIndex = { 0, 1, 2 };
303 for (
int i = 0; i < height; i++ )
312 feedback->
setProgress( 100.0 *
static_cast< double >( i ) / height );
318 scanLine->resetNoDataValue();
319 queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine->bits() );
321 memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i, 0 ), inputSize );
322 queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine->bits() );
324 memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i + 1, 0 ), inputSize );
325 queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine->bits() );
331 if ( i == inputBlock->height() - 1 )
333 scanLine->resetNoDataValue();
334 queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine->bits() );
338 queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, inputDataTypeSize * 1 , inputSize, inputBlock->bits( i + 1, 0 ) );
342 kernel( cl::EnqueueArgs( queue, cl::NDRange( width ) ), *scanLineBuffer[rowIndex[0]], *scanLineBuffer[rowIndex[1]], *scanLineBuffer[rowIndex[2]], resultLineBuffer, rasterParamsBuffer );
344 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, outputDataTypeSize * outputBlock->width(), outputBlock->bits( i, 0 ) );
345 std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
348 catch ( cl::Error &e )
360 double pixelValues[9] { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
361 bool isNoData[9] {
false,
false,
false,
false,
false,
false,
false,
false,
false };
363 for (
int row = 0; row < height; row++ )
365 for (
int col = 0; col < width; col++ )
373 else if ( row == height - 1 )
381 pixelValues[0] = inputBlock->valueAndNoData( iUp, 0, isNoData[0] );
382 pixelValues[1] = pixelValues[0];
383 isNoData[1] = isNoData[0];
384 pixelValues[2] = pixelValues[0];
385 isNoData[2] = isNoData[0];
387 pixelValues[3] = inputBlock->valueAndNoData( row, 0, isNoData[3] );
388 pixelValues[4] = pixelValues[3];
389 isNoData[4] = isNoData[3];
390 pixelValues[5] = pixelValues[3];
391 isNoData[5] = isNoData[3];
393 pixelValues[6] = inputBlock->valueAndNoData( iDown, 0, isNoData[6] );
394 pixelValues[7] = pixelValues[6];
395 isNoData[7] = isNoData[6];
396 pixelValues[8] = pixelValues[6];
397 isNoData[8] = isNoData[6];
402 pixelValues[0] = pixelValues[1];
403 pixelValues[1] = pixelValues[2];
404 pixelValues[3] = pixelValues[4];
405 pixelValues[4] = pixelValues[5];
406 pixelValues[6] = pixelValues[7];
407 pixelValues[7] = pixelValues[8];
408 isNoData[0] = isNoData[1];
409 isNoData[1] = isNoData[2];
410 isNoData[3] = isNoData[4];
411 isNoData[4] = isNoData[5];
412 isNoData[6] = isNoData[7];
413 isNoData[7] = isNoData[8];
417 if ( col < width - 1 )
419 pixelValues[2] = inputBlock->valueAndNoData( iUp, col + 1, isNoData[2] );
420 pixelValues[5] = inputBlock->valueAndNoData( row, col + 1, isNoData[5] );
421 pixelValues[8] = inputBlock->valueAndNoData( iDown, col + 1, isNoData[8] );
426 outputBlock->setColor( row, col, defaultNodataColor );
431 const double x22 = pixelValues[4];
433 const double x11 = isNoData[0] ? x22 : pixelValues[0];
434 const double x21 = isNoData[3] ? x22 : pixelValues[3];
435 const double x31 = isNoData[6] ? x22 : pixelValues[6];
436 const double x12 = isNoData[1] ? x22 : pixelValues[1];
438 const double x32 = isNoData[7] ? x22 : pixelValues[7];
439 const double x13 = isNoData[2] ? x22 : pixelValues[2];
440 const double x23 = isNoData[5] ? x22 : pixelValues[5];
441 const double x33 = isNoData[8] ? x22 : pixelValues[8];
444 const double derX = ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellXSize );
445 const double derY = ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellYSize );
450 if ( !mMultiDirectional )
454 = std::clamp( ( sin_altRadians_mul_254 - ( derY * cos_az_mul_cos_alt_mul_z_mul_254 - derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) / std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) ), 0.0, 255.0 );
460 const float xx = derX * derX;
461 const float yy = derY * derY;
462 const float xx_plus_yy = xx + yy;
464 if ( xx_plus_yy == 0.0 )
466 grayValue = std::clamp(
static_cast<float>( 1.0 + sin_altRadians_mul_254 ), 0.0f, 255.0f );
471 float val225_mul_127 = sin_altRadians_mul_127 + ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
472 val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
473 float val270_mul_127 = sin_altRadians_mul_127 - derX * cos_alt_mul_z_mul_127;
474 val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
475 float val315_mul_127 = sin_altRadians_mul_127 + ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
476 val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
477 float val360_mul_127 = sin_altRadians_mul_127 - derY * cos_alt_mul_z_mul_127;
478 val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
481 const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
482 const float weight_270 = xx;
483 const float weight_315 = xx_plus_yy - weight_225;
484 const float weight_360 = yy;
485 const float cang_mul_127 = ( ( weight_225 * val225_mul_127 + weight_270 * val270_mul_127 + weight_315 * val315_mul_127 + weight_360 * val360_mul_127 ) / xx_plus_yy )
486 / ( 1 + square_z * xx_plus_yy );
488 grayValue = std::clamp( 1.0f + cang_mul_127, 0.0f, 255.0f );
499 currentAlpha *= alphaBlock->value(
static_cast<qgssize>( row ) * width + col ) / 255.0;
504 outputBlock->setColor( row, col, qRgba( grayValue, grayValue, grayValue, 255 ) );
508 outputBlock->setColor( row, col, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
520 u
"%1 processing time for hillshade (%2 x %3 ): %4 ms"_s.arg( useOpenCL ? u
"OpenCL"_s : u
"CPU"_s )
523 .arg( std::chrono::duration_cast<std::chrono::milliseconds>( std::chrono::system_clock::now() - startTime ).count() ),
530 return outputBlock.release();
581 QDomNodeList elements = element.elementsByTagName( u
"sld:RasterSymbolizer"_s );
582 if ( elements.size() == 0 )
586 QDomElement rasterSymbolizerElem = elements.at( 0 ).toElement();
593 QDomElement channelSelectionElem = doc.createElement( u
"sld:ChannelSelection"_s );
594 elements = rasterSymbolizerElem.elementsByTagName( u
"sld:Opacity"_s );
595 if ( elements.size() != 0 )
597 rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
601 elements = rasterSymbolizerElem.elementsByTagName( u
"sld:Geometry"_s );
602 if ( elements.size() != 0 )
604 rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
608 rasterSymbolizerElem.insertBefore( channelSelectionElem, rasterSymbolizerElem.firstChild() );
613 QDomElement channelElem = doc.createElement( u
"sld:GrayChannel"_s );
614 channelSelectionElem.appendChild( channelElem );
617 QDomElement sourceChannelNameElem = doc.createElement( u
"sld:SourceChannelName"_s );
618 sourceChannelNameElem.appendChild( doc.createTextNode( QString::number( mBand ) ) );
619 channelElem.appendChild( sourceChannelNameElem );
623 QDomElement shadedReliefElem = doc.createElement( u
"sld:ShadedRelief"_s );
624 rasterSymbolizerElem.appendChild( shadedReliefElem );
627 QDomElement brightnessOnlyElem = doc.createElement( u
"sld:BrightnessOnly"_s );
628 brightnessOnlyElem.appendChild( doc.createTextNode( u
"true"_s ) );
629 shadedReliefElem.appendChild( brightnessOnlyElem );
632 QDomElement reliefFactorElem = doc.createElement( u
"sld:ReliefFactor"_s );
633 reliefFactorElem.appendChild( doc.createTextNode( QString::number(
zFactor() ) ) );
634 shadedReliefElem.appendChild( reliefFactorElem );
637 QDomElement altitudeVendorOptionElem = doc.createElement( u
"sld:VendorOption"_s );
638 altitudeVendorOptionElem.setAttribute( u
"name"_s, u
"altitude"_s );
639 altitudeVendorOptionElem.appendChild( doc.createTextNode( QString::number(
altitude() ) ) );
640 shadedReliefElem.appendChild( altitudeVendorOptionElem );
643 QDomElement azimutVendorOptionElem = doc.createElement( u
"sld:VendorOption"_s );
644 azimutVendorOptionElem.setAttribute( u
"name"_s, u
"azimuth"_s );
645 azimutVendorOptionElem.appendChild( doc.createTextNode( QString::number(
azimuth() ) ) );
646 shadedReliefElem.appendChild( azimutVendorOptionElem );
649 QDomElement multidirectionalVendorOptionElem = doc.createElement( u
"sld:VendorOption"_s );
650 multidirectionalVendorOptionElem.setAttribute( u
"name"_s, u
"multidirectional"_s );
651 multidirectionalVendorOptionElem.appendChild( doc.createTextNode( QString::number(
multiDirectional() ) ) );
652 shadedReliefElem.appendChild( multidirectionalVendorOptionElem );