QGIS API Documentation  3.13.0-Master (b73bd58cfb)
qgshillshaderenderer.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgshillshaderenderer.cpp
3  ---------------------------------
4  begin : May 2016
5  copyright : (C) 2016 by Nathan Woodrow
6  email : woodrow dot nathan at gmail dot com
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include <QColor>
19 
20 #include "qgshillshaderenderer.h"
21 #include "qgsrastertransparency.h"
22 #include "qgsrasterinterface.h"
23 #include "qgsrasterblock.h"
24 #include "qgsrectangle.h"
25 #include "qgsmessagelog.h"
26 #include <memory>
27 
28 #ifdef HAVE_OPENCL
29 #ifdef QGISDEBUG
30 #include <chrono>
31 #include "qgssettings.h"
32 #endif
33 #include "qgsexception.h"
34 #include "qgsopenclutils.h"
35 #endif
36 
37 QgsHillshadeRenderer::QgsHillshadeRenderer( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ):
38  QgsRasterRenderer( input, QStringLiteral( "hillshade" ) )
39  , mBand( band )
40  , mZFactor( 1 )
41  , mLightAngle( lightAngle )
42  , mLightAzimuth( lightAzimuth )
43  , mMultiDirectional( false )
44 {
45 
46 }
47 
49 {
50  QgsHillshadeRenderer *r = new QgsHillshadeRenderer( nullptr, mBand, mLightAzimuth, mLightAngle );
51  r->copyCommonProperties( this );
52 
53  r->setZFactor( mZFactor );
54  r->setMultiDirectional( mMultiDirectional );
55  return r;
56 }
57 
59 {
60  if ( elem.isNull() )
61  {
62  return nullptr;
63  }
64 
65  int band = elem.attribute( QStringLiteral( "band" ), QStringLiteral( "0" ) ).toInt();
66  double azimuth = elem.attribute( QStringLiteral( "azimuth" ), QStringLiteral( "315" ) ).toDouble();
67  double angle = elem.attribute( QStringLiteral( "angle" ), QStringLiteral( "45" ) ).toDouble();
68  double zFactor = elem.attribute( QStringLiteral( "zfactor" ), QStringLiteral( "1" ) ).toDouble();
69  bool multiDirectional = elem.attribute( QStringLiteral( "multidirection" ), QStringLiteral( "0" ) ).toInt();
70  QgsHillshadeRenderer *r = new QgsHillshadeRenderer( input, band, azimuth, angle );
71  r->readXml( elem );
72 
73  r->setZFactor( zFactor );
74  r->setMultiDirectional( multiDirectional );
75  return r;
76 }
77 
78 void QgsHillshadeRenderer::writeXml( QDomDocument &doc, QDomElement &parentElem ) const
79 {
80  if ( parentElem.isNull() )
81  {
82  return;
83  }
84 
85  QDomElement rasterRendererElem = doc.createElement( QStringLiteral( "rasterrenderer" ) );
86  _writeXml( doc, rasterRendererElem );
87 
88  rasterRendererElem.setAttribute( QStringLiteral( "band" ), mBand );
89  rasterRendererElem.setAttribute( QStringLiteral( "azimuth" ), QString::number( mLightAzimuth ) );
90  rasterRendererElem.setAttribute( QStringLiteral( "angle" ), QString::number( mLightAngle ) );
91  rasterRendererElem.setAttribute( QStringLiteral( "zfactor" ), QString::number( mZFactor ) );
92  rasterRendererElem.setAttribute( QStringLiteral( "multidirection" ), QString::number( mMultiDirectional ) );
93  parentElem.appendChild( rasterRendererElem );
94 }
95 
96 QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
97 {
98  Q_UNUSED( bandNo )
99  std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock() );
100  if ( !mInput )
101  {
102  QgsDebugMsg( QStringLiteral( "No input raster!" ) );
103  return outputBlock.release();
104  }
105 
106  std::shared_ptr< QgsRasterBlock > inputBlock( mInput->block( mBand, extent, width, height, feedback ) );
107 
108  if ( !inputBlock || inputBlock->isEmpty() )
109  {
110  QgsDebugMsg( QStringLiteral( "No raster data!" ) );
111  return outputBlock.release();
112  }
113 
114  std::shared_ptr< QgsRasterBlock > alphaBlock;
115 
116  if ( mAlphaBand > 0 && mBand != mAlphaBand )
117  {
118  alphaBlock.reset( mInput->block( mAlphaBand, extent, width, height, feedback ) );
119  if ( !alphaBlock || alphaBlock->isEmpty() )
120  {
121  // TODO: better to render without alpha
122  return outputBlock.release();
123  }
124  }
125  else if ( mAlphaBand > 0 )
126  {
127  alphaBlock = inputBlock;
128  }
129 
130  if ( !outputBlock->reset( Qgis::ARGB32_Premultiplied, width, height ) )
131  {
132  return outputBlock.release();
133  }
134 
135  // Starting the computation
136 
137  // Common pre-calculated values
138  float cellXSize = static_cast<float>( extent.width() ) / width;
139  float cellYSize = static_cast<float>( extent.height() ) / height;
140  float zenithRad = static_cast<float>( std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0 );
141  float azimuthRad = static_cast<float>( -1 * mLightAzimuth * M_PI / 180.0 );
142  float cosZenithRad = std::cos( zenithRad );
143  float sinZenithRad = std::sin( zenithRad );
144 
145  // For fast formula from GDAL DEM
146  float cos_alt_mul_z = cosZenithRad * static_cast<float>( mZFactor );
147  float cos_az_mul_cos_alt_mul_z = std::cos( azimuthRad ) * cos_alt_mul_z;
148  float sin_az_mul_cos_alt_mul_z = std::sin( azimuthRad ) * cos_alt_mul_z;
149  float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0f * cos_az_mul_cos_alt_mul_z;
150  float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0f * sin_az_mul_cos_alt_mul_z;
151  float square_z = static_cast<float>( mZFactor * mZFactor );
152  float sin_altRadians_mul_254 = 254.0f * sinZenithRad;
153 
154  // For multi directional
155  float sin_altRadians_mul_127 = 127.0f * sinZenithRad;
156  // 127.0 * std::cos(225.0 * M_PI / 180.0) = -32.87001872802012
157  float cos225_az_mul_cos_alt_mul_z_mul_127 = -32.87001872802012f * cos_alt_mul_z;
158  float cos_alt_mul_z_mul_127 = 127.0f * cos_alt_mul_z;
159 
160  const QRgb defaultNodataColor = renderColorForNodataPixel();
161 
162 #ifdef HAVE_OPENCL
163 
164  // Use OpenCL? For now OpenCL is enabled in the default configuration only
165  bool useOpenCL( QgsOpenClUtils::enabled()
168  && mAlphaBand <= 0
169  && inputBlock->dataTypeSize() <= 4 );
170  // Check for sources
171  QString source;
172  if ( useOpenCL )
173  {
174  source = QgsOpenClUtils::sourceFromBaseName( QStringLiteral( "hillshade_renderer" ) );
175  if ( source.isEmpty() )
176  {
177  useOpenCL = false;
178  QgsMessageLog::logMessage( QObject::tr( "Error loading OpenCL program source from path %1" ).arg( QgsOpenClUtils::sourcePath() ), QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
179  }
180  }
181 
182 #ifdef QGISDEBUG
183  std::chrono::time_point<std::chrono::system_clock> startTime( std::chrono::system_clock::now() );
184 #endif
185 
186  if ( useOpenCL )
187  {
188 
189  try
190  {
191  std::size_t inputDataTypeSize = inputBlock->dataTypeSize();
192  std::size_t outputDataTypeSize = outputBlock->dataTypeSize();
193  // Buffer scanline, 1px height, 2px wider
194  QString typeName;
195  switch ( inputBlock->dataType() )
196  {
197  case Qgis::DataType::Byte:
198  typeName = QStringLiteral( "unsigned char" );
199  break;
200  case Qgis::DataType::UInt16:
201  typeName = QStringLiteral( "unsigned int" );
202  break;
203  case Qgis::DataType::Int16:
204  typeName = QStringLiteral( "short" );
205  break;
206  case Qgis::DataType::UInt32:
207  typeName = QStringLiteral( "unsigned int" );
208  break;
209  case Qgis::DataType::Int32:
210  typeName = QStringLiteral( "int" );
211  break;
212  case Qgis::DataType::Float32:
213  typeName = QStringLiteral( "float" );
214  break;
215  default:
216  throw QgsException( QStringLiteral( "Unsupported data type for OpenCL processing." ) );
217  }
218 
219  if ( inputBlock->dataType() != Qgis::DataType::Float32 )
220  {
221  source.replace( QStringLiteral( "__global float *scanLine" ), QStringLiteral( "__global %1 *scanLine" ).arg( typeName ) );
222  }
223 
224  // Data type for input is Float32 (4 bytes)
225  std::size_t scanLineWidth( inputBlock->width() + 2 );
226  std::size_t inputSize( inputDataTypeSize * inputBlock->width() );
227 
228  // CL buffers are also 2px wider for nodata initial and final columns
229  std::size_t bufferWidth( width + 2 );
230  std::size_t bufferSize( inputDataTypeSize * bufferWidth );
231 
232  // Buffer scanlines, 1px height, 2px wider
233  // Data type for input is Float32 (4 bytes)
234  // keep only three scanlines in memory at a time, make room for initial and final nodata
235  std::unique_ptr<QgsRasterBlock> scanLine = qgis::make_unique<QgsRasterBlock>( inputBlock->dataType(), scanLineWidth, 1 );
236  // Note: output block is not 2px wider and it is an image
237  // Prepare context and queue
238  cl::Context ctx = QgsOpenClUtils::context();
239  cl::CommandQueue queue = QgsOpenClUtils::commandQueue();
240 
241  // Cast to float (because double just crashes on some GPUs)
242  std::vector<float> rasterParams;
243  rasterParams.push_back( inputBlock->noDataValue() );
244  rasterParams.push_back( outputBlock->noDataValue() );
245  rasterParams.push_back( mZFactor );
246  rasterParams.push_back( cellXSize );
247  rasterParams.push_back( cellYSize );
248  rasterParams.push_back( static_cast<float>( mOpacity ) ); // 5
249 
250  // For fast formula from GDAL DEM
251  rasterParams.push_back( cos_az_mul_cos_alt_mul_z_mul_254 ); // 6
252  rasterParams.push_back( sin_az_mul_cos_alt_mul_z_mul_254 ); // 7
253  rasterParams.push_back( square_z ); // 8
254  rasterParams.push_back( sin_altRadians_mul_254 ); // 9
255 
256  // For multidirectional fast formula
257  rasterParams.push_back( sin_altRadians_mul_127 ); // 10
258  rasterParams.push_back( cos225_az_mul_cos_alt_mul_z_mul_127 ); // 11
259  rasterParams.push_back( cos_alt_mul_z_mul_127 ); // 12
260 
261  // Default color for nodata (BGR components)
262  rasterParams.push_back( static_cast<float>( qBlue( defaultNodataColor ) ) ); // 13
263  rasterParams.push_back( static_cast<float>( qGreen( defaultNodataColor ) ) ); // 14
264  rasterParams.push_back( static_cast<float>( qRed( defaultNodataColor ) ) ); // 15
265  rasterParams.push_back( static_cast<float>( qAlpha( defaultNodataColor ) ) / 255.0f ); // 16
266 
267  // Whether use multidirectional
268  rasterParams.push_back( static_cast<float>( mMultiDirectional ) ); // 17
269 
270  cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr );
271  cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
272  cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
273  cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
274  cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer};
275  // Note that result buffer is an image
276  cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, outputDataTypeSize * width, nullptr, nullptr );
277 
278  static std::map<Qgis::DataType, cl::Program> programCache;
279  cl::Program program = programCache[inputBlock->dataType()];
280  if ( ! program.get() )
281  {
282  // Create a program from the kernel source
283  programCache[inputBlock->dataType()] = QgsOpenClUtils::buildProgram( source, QgsOpenClUtils::ExceptionBehavior::Throw );
284  program = programCache[inputBlock->dataType()];
285  }
286 
287  // Disable program cache when developing and testing cl program
288  // program = QgsOpenClUtils::buildProgram( ctx, source, QgsOpenClUtils::ExceptionBehavior::Throw );
289 
290  // Create the OpenCL kernel
291  auto kernel = cl::KernelFunctor <
292  cl::Buffer &,
293  cl::Buffer &,
294  cl::Buffer &,
295  cl::Buffer &,
296  cl::Buffer &
297  > ( program, "processNineCellWindow" );
298 
299 
300  // Rotating buffer index
301  std::vector<int> rowIndex = {0, 1, 2};
302 
303  for ( int i = 0; i < height; i++ )
304  {
305  if ( feedback && feedback->isCanceled() )
306  {
307  break;
308  }
309 
310  if ( feedback )
311  {
312  feedback->setProgress( 100.0 * static_cast< double >( i ) / height );
313  }
314 
315  if ( i == 0 )
316  {
317  // Fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
318  scanLine->resetNoDataValue();
319  queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) );
320  // Read first row
321  memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i, 0 ), inputSize );
322  queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); // row 0
323  // Second row
324  memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i + 1, 0 ), inputSize );
325  queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); //
326  }
327  else
328  {
329  // Normally fetch only scanLine3 and move forward one row
330  // Read scanline 3, fill the last row with nodata values if it's the last iteration
331  if ( i == inputBlock->height() - 1 )
332  {
333  scanLine->resetNoDataValue();
334  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine->bits( ) );
335  }
336  else // Overwrite from input, skip first and last
337  {
338  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, inputDataTypeSize * 1 /* offset 1 */, inputSize, inputBlock->bits( i + 1, 0 ) );
339  }
340  }
341 
342  kernel( cl::EnqueueArgs(
343  queue,
344  cl::NDRange( width )
345  ),
346  *scanLineBuffer[rowIndex[0]],
347  *scanLineBuffer[rowIndex[1]],
348  *scanLineBuffer[rowIndex[2]],
349  resultLineBuffer,
350  rasterParamsBuffer
351  );
352 
353  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, outputDataTypeSize * outputBlock->width( ), outputBlock->bits( i, 0 ) );
354  std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
355  }
356  }
357  catch ( cl::Error &e )
358  {
359  QgsMessageLog::logMessage( QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what( ) ).arg( QgsOpenClUtils::errorText( e.err( ) ) ),
362  QgsMessageLog::logMessage( QObject::tr( "OpenCL has been disabled, you can re-enable it in the options dialog." ),
364  }
365 
366  } // End of OpenCL processing path
367  else // Use the CPU and the original algorithm
368  {
369 
370 #endif
371 
372  for ( qgssize i = 0; i < static_cast<qgssize>( height ); i++ )
373  {
374 
375  for ( qgssize j = 0; j < static_cast<qgssize>( width ); j++ )
376  {
377 
378  if ( inputBlock->isNoData( i, j ) )
379  {
380  outputBlock->setColor( static_cast<int>( i ), static_cast<int>( j ), defaultNodataColor );
381  continue;
382  }
383 
384  qgssize iUp, iDown, jLeft, jRight;
385  if ( i == 0 )
386  {
387  iUp = i;
388  iDown = i + 1;
389  }
390  else if ( i < static_cast<qgssize>( height ) - 1 )
391  {
392  iUp = i - 1;
393  iDown = i + 1;
394  }
395  else
396  {
397  iUp = i - 1;
398  iDown = i;
399  }
400 
401  if ( j == 0 )
402  {
403  jLeft = j;
404  jRight = j + 1;
405  }
406  else if ( j < static_cast<qgssize>( width ) - 1 )
407  {
408  jLeft = j - 1;
409  jRight = j + 1;
410  }
411  else
412  {
413  jLeft = j - 1;
414  jRight = j;
415  }
416 
417  double x11;
418  double x21;
419  double x31;
420  double x12;
421  double x22; // Working cell
422  double x32;
423  double x13;
424  double x23;
425  double x33;
426 
427  // This is center cell. It is not nodata. Use this in place of nodata neighbors
428  x22 = inputBlock->value( i, j );
429 
430  x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft );
431  x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft );
432  x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft );
433 
434  x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j );
435  // x22
436  x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j );
437 
438  x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight );
439  x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight );
440  x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight );
441 
442  double derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize );
443  double derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize );
444 
445  // Fast formula
446 
447  double grayValue;
448  if ( !mMultiDirectional )
449  {
450  // Standard single direction hillshade
451  grayValue = qBound( 0.0, ( sin_altRadians_mul_254 -
452  ( derY * cos_az_mul_cos_alt_mul_z_mul_254 -
453  derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) /
454  std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) )
455  , 255.0 );
456  }
457  else
458  {
459  // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
460  // Fast formula from GDAL DEM
461  const float xx = derX * derX;
462  const float yy = derY * derY;
463  const float xx_plus_yy = xx + yy;
464  // Flat?
465  if ( xx_plus_yy == 0.0 )
466  {
467  grayValue = qBound( 0.0f, static_cast<float>( 1.0 + sin_altRadians_mul_254 ), 255.0f );
468  }
469  else
470  {
471  // ... then the shade value from different azimuth
472  float val225_mul_127 = sin_altRadians_mul_127 +
473  ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
474  val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
475  float val270_mul_127 = sin_altRadians_mul_127 -
476  derX * cos_alt_mul_z_mul_127;
477  val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
478  float val315_mul_127 = sin_altRadians_mul_127 +
479  ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
480  val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
481  float val360_mul_127 = sin_altRadians_mul_127 -
482  derY * cos_alt_mul_z_mul_127;
483  val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
484 
485  // ... then the weighted shading
486  const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
487  const float weight_270 = xx;
488  const float weight_315 = xx_plus_yy - weight_225;
489  const float weight_360 = yy;
490  const float cang_mul_127 = (
491  ( weight_225 * val225_mul_127 +
492  weight_270 * val270_mul_127 +
493  weight_315 * val315_mul_127 +
494  weight_360 * val360_mul_127 ) / xx_plus_yy ) /
495  ( 1 + square_z * xx_plus_yy );
496 
497  grayValue = qBound( 0.0f, 1.0f + cang_mul_127, 255.0f );
498  }
499  }
500 
501  double currentAlpha = mOpacity;
502  if ( mRasterTransparency )
503  {
504  currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0;
505  }
506  if ( mAlphaBand > 0 )
507  {
508  currentAlpha *= alphaBlock->value( i ) / 255.0;
509  }
510 
511  if ( qgsDoubleNear( currentAlpha, 1.0 ) )
512  {
513  outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) );
514  }
515  else
516  {
517  outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
518  }
519  }
520  }
521 
522 #ifdef HAVE_OPENCL
523  } // End of switch in case OpenCL is not available or enabled
524 
525 #ifdef QGISDEBUG
526  if ( QgsSettings().value( QStringLiteral( "Map/logCanvasRefreshEvent" ), false ).toBool() )
527  {
528  QgsMessageLog::logMessage( QStringLiteral( "%1 processing time for hillshade (%2 x %3 ): %4 ms" )
529  .arg( useOpenCL ? QStringLiteral( "OpenCL" ) : QStringLiteral( "CPU" ) )
530  .arg( width )
531  .arg( height )
532  .arg( std::chrono::duration_cast<std::chrono::milliseconds>( std::chrono::system_clock::now() - startTime ).count() ),
533  tr( "Rendering" ) );
534  }
535 #endif
536 #endif
537 
538  return outputBlock.release();
539 }
540 
542 {
543  QList<int> bandList;
544  if ( mBand != -1 )
545  {
546  bandList << mBand;
547  }
548  return bandList;
549 
550 }
551 
553 {
554  if ( bandNo > mInput->bandCount() || bandNo <= 0 )
555  {
556  return;
557  }
558  mBand = bandNo;
559 }
560 
561 double QgsHillshadeRenderer::calcFirstDerX( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
562 {
563  Q_UNUSED( x12 )
564  Q_UNUSED( x22 )
565  Q_UNUSED( x32 )
566  return ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellsize );
567 }
568 
569 double QgsHillshadeRenderer::calcFirstDerY( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
570 {
571  Q_UNUSED( x21 )
572  Q_UNUSED( x22 )
573  Q_UNUSED( x23 )
574  return ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellsize );
575 }
576 
577 void QgsHillshadeRenderer::toSld( QDomDocument &doc, QDomElement &element, const QgsStringMap &props ) const
578 {
579  // create base structure
580  QgsRasterRenderer::toSld( doc, element, props );
581 
582  // look for RasterSymbolizer tag
583  QDomNodeList elements = element.elementsByTagName( QStringLiteral( "sld:RasterSymbolizer" ) );
584  if ( elements.size() == 0 )
585  return;
586 
587  // there SHOULD be only one
588  QDomElement rasterSymbolizerElem = elements.at( 0 ).toElement();
589 
590  // add Channel Selection tags (if band is not default 1)
591  // Need to insert channelSelection in the correct sequence as in SLD standard e.g.
592  // after opacity or geometry or as first element after sld:RasterSymbolizer
593  if ( band() != 1 )
594  {
595  QDomElement channelSelectionElem = doc.createElement( QStringLiteral( "sld:ChannelSelection" ) );
596  elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral( "sld:Opacity" ) );
597  if ( elements.size() != 0 )
598  {
599  rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
600  }
601  else
602  {
603  elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral( "sld:Geometry" ) );
604  if ( elements.size() != 0 )
605  {
606  rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
607  }
608  else
609  {
610  rasterSymbolizerElem.insertBefore( channelSelectionElem, rasterSymbolizerElem.firstChild() );
611  }
612  }
613 
614  // for gray band
615  QDomElement channelElem = doc.createElement( QStringLiteral( "sld:GrayChannel" ) );
616  channelSelectionElem.appendChild( channelElem );
617 
618  // set band
619  QDomElement sourceChannelNameElem = doc.createElement( QStringLiteral( "sld:SourceChannelName" ) );
620  sourceChannelNameElem.appendChild( doc.createTextNode( QString::number( band() ) ) );
621  channelElem.appendChild( sourceChannelNameElem );
622  }
623 
624  // add ShadedRelief tag
625  QDomElement shadedReliefElem = doc.createElement( QStringLiteral( "sld:ShadedRelief" ) );
626  rasterSymbolizerElem.appendChild( shadedReliefElem );
627 
628  // brightnessOnly tag
629  QDomElement brightnessOnlyElem = doc.createElement( QStringLiteral( "sld:BrightnessOnly" ) );
630  brightnessOnlyElem.appendChild( doc.createTextNode( QStringLiteral( "true" ) ) );
631  shadedReliefElem.appendChild( brightnessOnlyElem );
632 
633  // ReliefFactor tag
634  QDomElement reliefFactorElem = doc.createElement( QStringLiteral( "sld:ReliefFactor" ) );
635  reliefFactorElem.appendChild( doc.createTextNode( QString::number( zFactor() ) ) );
636  shadedReliefElem.appendChild( reliefFactorElem );
637 
638  // altitude VendorOption tag
639  QDomElement altitudeVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
640  altitudeVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "altitude" ) );
641  altitudeVendorOptionElem.appendChild( doc.createTextNode( QString::number( altitude() ) ) );
642  shadedReliefElem.appendChild( altitudeVendorOptionElem );
643 
644  // azimuth VendorOption tag
645  QDomElement azimutVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
646  azimutVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "azimuth" ) );
647  azimutVendorOptionElem.appendChild( doc.createTextNode( QString::number( azimuth() ) ) );
648  shadedReliefElem.appendChild( azimutVendorOptionElem );
649 
650  // multidirectional VendorOption tag
651  QDomElement multidirectionalVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
652  multidirectionalVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "multidirectional" ) );
653  multidirectionalVendorOptionElem.appendChild( doc.createTextNode( QString::number( multiDirectional() ) ) );
654  shadedReliefElem.appendChild( multidirectionalVendorOptionElem );
655 }
virtual int bandCount() const =0
Gets number of bands.
static QString sourcePath()
Returns the base path to OpenCL program directory.
A rectangle specified with double values.
Definition: qgsrectangle.h:41
bool isEmpty() const
True if there are no entries in the pixel lists except the nodata value.
static void setEnabled(bool enabled)
Set the OpenCL user setting to enabled.
static QString sourceFromBaseName(const QString &baseName)
Returns the full path to a an OpenCL source file from the baseName (&#39;.cl&#39; extension is automatically ...
virtual void toSld(QDomDocument &doc, QDomElement &element, const QgsStringMap &props=QgsStringMap()) const
Used from subclasses to create SLD Rule elements following SLD v1.0 specs.
void setMultiDirectional(bool isMultiDirectional)
Sets whether to render using a multi-directional hillshade algorithm.
This class is a composition of two QSettings instances:
Definition: qgssettings.h:61
virtual QgsRectangle extent() const
Gets the extent of the interface.
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:315
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:64
virtual QgsRasterInterface * input() const
Current input.
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
QMap< QString, QString > QgsStringMap
Definition: qgis.h:694
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:786
static cl::Context context()
Context factory.
void setBand(int bandNo)
Sets the band used by the renderer.
QgsRasterTransparency * mRasterTransparency
Raster transparency per color or value. Overwrites global alpha value.
void setZFactor(double zfactor)
Set the Z scaling factor of the result image.
int band() const
Returns the band used by the renderer.
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:116
QList< int > usesBands() const override
Returns a list of band numbers used by the renderer.
Raster data container.
bool multiDirectional() const
Returns true if the renderer is using multi-directional hillshading.
virtual QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr)=0
Read block of data using given extent and size.
const QString & typeName
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
void _writeXml(QDomDocument &doc, QDomElement &rasterRendererElem) const
Write upper class info into rasterrenderer element (called by writeXml method of subclasses) ...
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
void copyCommonProperties(const QgsRasterRenderer *other, bool copyMinMaxOrigin=true)
Copies common properties like opacity / transparency data from other renderer.
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
A renderer for generating live hillshade models.
QgsHillshadeRenderer(QgsRasterInterface *input, int band, double lightAzimuth, double lightAltitude)
A renderer for generating live hillshade models.
int mAlphaBand
Read alpha value from band.
void readXml(const QDomElement &rendererElem) override
Sets base class members from xml. Usually called from create() methods of subclasses.
Base class for processing filters like renderers, reprojector, resampler etc.
int alphaValue(double value, int globalTransparency=255) const
Returns the transparency value for a single value pixel.
static QString errorText(const int errorCode)
Returns a string representation from an OpenCL errorCode.
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition: qgis.h:703
QgsHillshadeRenderer * clone() const override
Clone itself, create deep copy.
static Q_DECL_DEPRECATED cl::Program buildProgram(const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior=Catch)
Build the program from source in the given context and depending on exceptionBehavior can throw or ca...
void toSld(QDomDocument &doc, QDomElement &element, const QgsStringMap &props=QgsStringMap()) const override
Used from subclasses to create SLD Rule elements following SLD v1.0 specs.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:55
static QLatin1String LOGMESSAGE_TAG
OpenCL string for message logs.
QRgb renderColorForNodataPixel() const
Returns the color for the renderer to use to represent nodata pixels.
double altitude() const
Returns the angle of the light source over the raster.
void writeXml(QDomDocument &doc, QDomElement &parentElem) const override
Write base class members to xml.
double mOpacity
Global alpha value (0-1)
double zFactor() const
Returns the Z scaling factor.
double azimuth() const
Returns the direction of the light over the raster between 0-360.
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
static QgsRasterRenderer * create(const QDomElement &elem, QgsRasterInterface *input)
Factory method to create a new renderer.
QgsRasterInterface * mInput
Feedback object tailored for raster block reading.
Defines a QGIS exception class.
Definition: qgsexception.h:34
Raster renderer pipe that applies colors to a raster.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:209