QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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  QRgb defaultNodataColor = NODATA_COLOR;
161 
162 
163 #ifdef HAVE_OPENCL
164 
165  // Use OpenCL? For now OpenCL is enabled in the default configuration only
166  bool useOpenCL( QgsOpenClUtils::enabled()
169  && mAlphaBand <= 0
170  && inputBlock->dataTypeSize() <= 4 );
171  // Check for sources
172  QString source;
173  if ( useOpenCL )
174  {
175  source = QgsOpenClUtils::sourceFromBaseName( QStringLiteral( "hillshade_renderer" ) );
176  if ( source.isEmpty() )
177  {
178  useOpenCL = false;
179  QgsMessageLog::logMessage( QObject::tr( "Error loading OpenCL program source from path %1" ).arg( QgsOpenClUtils::sourcePath() ), QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
180  }
181  }
182 
183 #ifdef QGISDEBUG
184  std::chrono::time_point<std::chrono::system_clock> startTime( std::chrono::system_clock::now() );
185 #endif
186 
187  if ( useOpenCL )
188  {
189 
190  try
191  {
192  std::size_t inputDataTypeSize = inputBlock->dataTypeSize();
193  std::size_t outputDataTypeSize = outputBlock->dataTypeSize();
194  // Buffer scanline, 1px height, 2px wider
195  QString typeName;
196  switch ( inputBlock->dataType() )
197  {
198  case Qgis::DataType::Byte:
199  typeName = QStringLiteral( "unsigned char" );
200  break;
201  case Qgis::DataType::UInt16:
202  typeName = QStringLiteral( "unsigned int" );
203  break;
204  case Qgis::DataType::Int16:
205  typeName = QStringLiteral( "short" );
206  break;
207  case Qgis::DataType::UInt32:
208  typeName = QStringLiteral( "unsigned int" );
209  break;
210  case Qgis::DataType::Int32:
211  typeName = QStringLiteral( "int" );
212  break;
213  case Qgis::DataType::Float32:
214  typeName = QStringLiteral( "float" );
215  break;
216  default:
217  throw QgsException( QStringLiteral( "Unsupported data type for OpenCL processing." ) );
218  }
219 
220  if ( inputBlock->dataType() != Qgis::DataType::Float32 )
221  {
222  source.replace( QStringLiteral( "__global float *scanLine" ), QStringLiteral( "__global %1 *scanLine" ).arg( typeName ) );
223  }
224 
225  // Data type for input is Float32 (4 bytes)
226  std::size_t scanLineWidth( inputBlock->width() + 2 );
227  std::size_t inputSize( inputDataTypeSize * inputBlock->width() );
228 
229  // CL buffers are also 2px wider for nodata initial and final columns
230  std::size_t bufferWidth( width + 2 );
231  std::size_t bufferSize( inputDataTypeSize * bufferWidth );
232 
233  // Buffer scanlines, 1px height, 2px wider
234  // Data type for input is Float32 (4 bytes)
235  // keep only three scanlines in memory at a time, make room for initial and final nodata
236  std::unique_ptr<QgsRasterBlock> scanLine = qgis::make_unique<QgsRasterBlock>( inputBlock->dataType(), scanLineWidth, 1 );
237  // Note: output block is not 2px wider and it is an image
238  // Prepare context and queue
239  cl::Context ctx = QgsOpenClUtils::context();
240  cl::CommandQueue queue = QgsOpenClUtils::commandQueue();
241 
242  // Cast to float (because double just crashes on some GPUs)
243  std::vector<float> rasterParams;
244  rasterParams.push_back( inputBlock->noDataValue() );
245  rasterParams.push_back( outputBlock->noDataValue() );
246  rasterParams.push_back( mZFactor );
247  rasterParams.push_back( cellXSize );
248  rasterParams.push_back( cellYSize );
249  rasterParams.push_back( static_cast<float>( mOpacity ) ); // 5
250 
251  // For fast formula from GDAL DEM
252  rasterParams.push_back( cos_az_mul_cos_alt_mul_z_mul_254 ); // 6
253  rasterParams.push_back( sin_az_mul_cos_alt_mul_z_mul_254 ); // 7
254  rasterParams.push_back( square_z ); // 8
255  rasterParams.push_back( sin_altRadians_mul_254 ); // 9
256 
257  // For multidirectional fast formula
258  rasterParams.push_back( sin_altRadians_mul_127 ); // 10
259  rasterParams.push_back( cos225_az_mul_cos_alt_mul_z_mul_127 ); // 11
260  rasterParams.push_back( cos_alt_mul_z_mul_127 ); // 12
261 
262  // Default color for nodata (BGR components)
263  rasterParams.push_back( static_cast<float>( qBlue( defaultNodataColor ) ) ); // 13
264  rasterParams.push_back( static_cast<float>( qGreen( defaultNodataColor ) ) ); // 14
265  rasterParams.push_back( static_cast<float>( qRed( defaultNodataColor ) ) ); // 15
266  rasterParams.push_back( static_cast<float>( qAlpha( defaultNodataColor ) ) / 255.0f ); // 16
267 
268  // Whether use multidirectional
269  rasterParams.push_back( static_cast<float>( mMultiDirectional ) ); // 17
270 
271  cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr );
272  cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
273  cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
274  cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
275  cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer};
276  // Note that result buffer is an image
277  cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, outputDataTypeSize * width, nullptr, nullptr );
278 
279  static std::map<Qgis::DataType, cl::Program> programCache;
280  cl::Program program = programCache[inputBlock->dataType()];
281  if ( ! program.get() )
282  {
283  // Create a program from the kernel source
284  programCache[inputBlock->dataType()] = QgsOpenClUtils::buildProgram( source, QgsOpenClUtils::ExceptionBehavior::Throw );
285  program = programCache[inputBlock->dataType()];
286  }
287 
288  // Disable program cache when developing and testing cl program
289  // program = QgsOpenClUtils::buildProgram( ctx, source, QgsOpenClUtils::ExceptionBehavior::Throw );
290 
291  // Create the OpenCL kernel
292  auto kernel = cl::KernelFunctor <
293  cl::Buffer &,
294  cl::Buffer &,
295  cl::Buffer &,
296  cl::Buffer &,
297  cl::Buffer &
298  > ( program, "processNineCellWindow" );
299 
300 
301  // Rotating buffer index
302  std::vector<int> rowIndex = {0, 1, 2};
303 
304  for ( int i = 0; i < height; i++ )
305  {
306  if ( feedback && feedback->isCanceled() )
307  {
308  break;
309  }
310 
311  if ( feedback )
312  {
313  feedback->setProgress( 100.0 * static_cast< double >( i ) / height );
314  }
315 
316  if ( i == 0 )
317  {
318  // Fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
319  scanLine->resetNoDataValue();
320  queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) );
321  // Read first row
322  memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i, 0 ), inputSize );
323  queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); // row 0
324  // Second row
325  memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i + 1, 0 ), inputSize );
326  queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); //
327  }
328  else
329  {
330  // Normally fetch only scanLine3 and move forward one row
331  // Read scanline 3, fill the last row with nodata values if it's the last iteration
332  if ( i == inputBlock->height() - 1 )
333  {
334  scanLine->resetNoDataValue();
335  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine->bits( ) );
336  }
337  else // Overwrite from input, skip first and last
338  {
339  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, inputDataTypeSize * 1 /* offset 1 */, inputSize, inputBlock->bits( i + 1, 0 ) );
340  }
341  }
342 
343  kernel( cl::EnqueueArgs(
344  queue,
345  cl::NDRange( width )
346  ),
347  *scanLineBuffer[rowIndex[0]],
348  *scanLineBuffer[rowIndex[1]],
349  *scanLineBuffer[rowIndex[2]],
350  resultLineBuffer,
351  rasterParamsBuffer
352  );
353 
354  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, outputDataTypeSize * outputBlock->width( ), outputBlock->bits( i, 0 ) );
355  std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
356  }
357  }
358  catch ( cl::Error &e )
359  {
360  QgsMessageLog::logMessage( QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what( ) ).arg( QgsOpenClUtils::errorText( e.err( ) ) ),
363  QgsMessageLog::logMessage( QObject::tr( "OpenCL has been disabled, you can re-enable it in the options dialog." ),
365  }
366 
367  } // End of OpenCL processing path
368  else // Use the CPU and the original algorithm
369  {
370 
371 #endif
372 
373  for ( qgssize i = 0; i < static_cast<qgssize>( height ); i++ )
374  {
375 
376  for ( qgssize j = 0; j < static_cast<qgssize>( width ); j++ )
377  {
378 
379  if ( inputBlock->isNoData( i, j ) )
380  {
381  outputBlock->setColor( static_cast<int>( i ), static_cast<int>( j ), defaultNodataColor );
382  continue;
383  }
384 
385  qgssize iUp, iDown, jLeft, jRight;
386  if ( i == 0 )
387  {
388  iUp = i;
389  iDown = i + 1;
390  }
391  else if ( i < static_cast<qgssize>( height ) - 1 )
392  {
393  iUp = i - 1;
394  iDown = i + 1;
395  }
396  else
397  {
398  iUp = i - 1;
399  iDown = i;
400  }
401 
402  if ( j == 0 )
403  {
404  jLeft = j;
405  jRight = j + 1;
406  }
407  else if ( j < static_cast<qgssize>( width ) - 1 )
408  {
409  jLeft = j - 1;
410  jRight = j + 1;
411  }
412  else
413  {
414  jLeft = j - 1;
415  jRight = j;
416  }
417 
418  double x11;
419  double x21;
420  double x31;
421  double x12;
422  double x22; // Working cell
423  double x32;
424  double x13;
425  double x23;
426  double x33;
427 
428  // This is center cell. It is not nodata. Use this in place of nodata neighbors
429  x22 = inputBlock->value( i, j );
430 
431  x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft );
432  x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft );
433  x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft );
434 
435  x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j );
436  // x22
437  x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j );
438 
439  x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight );
440  x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight );
441  x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight );
442 
443  double derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize );
444  double derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize );
445 
446  // Fast formula
447 
448  double grayValue;
449  if ( !mMultiDirectional )
450  {
451  // Standard single direction hillshade
452  grayValue = qBound( 0.0, ( sin_altRadians_mul_254 -
453  ( derY * cos_az_mul_cos_alt_mul_z_mul_254 -
454  derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) /
455  std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) )
456  , 255.0 );
457  }
458  else
459  {
460  // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
461  // Fast formula from GDAL DEM
462  const float xx = derX * derX;
463  const float yy = derY * derY;
464  const float xx_plus_yy = xx + yy;
465  // Flat?
466  if ( xx_plus_yy == 0.0 )
467  {
468  grayValue = qBound( 0.0f, static_cast<float>( 1.0 + sin_altRadians_mul_254 ), 255.0f );
469  }
470  else
471  {
472  // ... then the shade value from different azimuth
473  float val225_mul_127 = sin_altRadians_mul_127 +
474  ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
475  val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
476  float val270_mul_127 = sin_altRadians_mul_127 -
477  derX * cos_alt_mul_z_mul_127;
478  val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
479  float val315_mul_127 = sin_altRadians_mul_127 +
480  ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
481  val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
482  float val360_mul_127 = sin_altRadians_mul_127 -
483  derY * cos_alt_mul_z_mul_127;
484  val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
485 
486  // ... then the weighted shading
487  const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
488  const float weight_270 = xx;
489  const float weight_315 = xx_plus_yy - weight_225;
490  const float weight_360 = yy;
491  const float cang_mul_127 = (
492  ( weight_225 * val225_mul_127 +
493  weight_270 * val270_mul_127 +
494  weight_315 * val315_mul_127 +
495  weight_360 * val360_mul_127 ) / xx_plus_yy ) /
496  ( 1 + square_z * xx_plus_yy );
497 
498  grayValue = qBound( 0.0f, 1.0f + cang_mul_127, 255.0f );
499  }
500  }
501 
502  double currentAlpha = mOpacity;
503  if ( mRasterTransparency )
504  {
505  currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0;
506  }
507  if ( mAlphaBand > 0 )
508  {
509  currentAlpha *= alphaBlock->value( i ) / 255.0;
510  }
511 
512  if ( qgsDoubleNear( currentAlpha, 1.0 ) )
513  {
514  outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) );
515  }
516  else
517  {
518  outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
519  }
520  }
521  }
522 
523 #ifdef HAVE_OPENCL
524  } // End of switch in case OpenCL is not available or enabled
525 
526 #ifdef QGISDEBUG
527  if ( QgsSettings().value( QStringLiteral( "Map/logCanvasRefreshEvent" ), false ).toBool() )
528  {
529  QgsMessageLog::logMessage( QStringLiteral( "%1 processing time for hillshade (%2 x %3 ): %4 ms" )
530  .arg( useOpenCL ? QStringLiteral( "OpenCL" ) : QStringLiteral( "CPU" ) )
531  .arg( width )
532  .arg( height )
533  .arg( std::chrono::duration_cast<std::chrono::milliseconds>( std::chrono::system_clock::now() - startTime ).count() ),
534  tr( "Rendering" ) );
535  }
536 #endif
537 #endif
538 
539  return outputBlock.release();
540 }
541 
543 {
544  QList<int> bandList;
545  if ( mBand != -1 )
546  {
547  bandList << mBand;
548  }
549  return bandList;
550 
551 }
552 
554 {
555  if ( bandNo > mInput->bandCount() || bandNo <= 0 )
556  {
557  return;
558  }
559  mBand = bandNo;
560 }
561 
562 double QgsHillshadeRenderer::calcFirstDerX( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
563 {
564  Q_UNUSED( x12 )
565  Q_UNUSED( x22 )
566  Q_UNUSED( x32 )
567  return ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellsize );
568 }
569 
570 double QgsHillshadeRenderer::calcFirstDerY( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
571 {
572  Q_UNUSED( x21 )
573  Q_UNUSED( x22 )
574  Q_UNUSED( x23 )
575  return ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellsize );
576 }
577 
578 void QgsHillshadeRenderer::toSld( QDomDocument &doc, QDomElement &element, const QgsStringMap &props ) const
579 {
580  // create base structure
581  QgsRasterRenderer::toSld( doc, element, props );
582 
583  // look for RasterSymbolizer tag
584  QDomNodeList elements = element.elementsByTagName( QStringLiteral( "sld:RasterSymbolizer" ) );
585  if ( elements.size() == 0 )
586  return;
587 
588  // there SHOULD be only one
589  QDomElement rasterSymbolizerElem = elements.at( 0 ).toElement();
590 
591  // add Channel Selection tags (if band is not default 1)
592  // Need to insert channelSelection in the correct sequence as in SLD standard e.g.
593  // after opacity or geometry or as first element after sld:RasterSymbolizer
594  if ( band() != 1 )
595  {
596  QDomElement channelSelectionElem = doc.createElement( QStringLiteral( "sld:ChannelSelection" ) );
597  elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral( "sld:Opacity" ) );
598  if ( elements.size() != 0 )
599  {
600  rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
601  }
602  else
603  {
604  elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral( "sld:Geometry" ) );
605  if ( elements.size() != 0 )
606  {
607  rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
608  }
609  else
610  {
611  rasterSymbolizerElem.insertBefore( channelSelectionElem, rasterSymbolizerElem.firstChild() );
612  }
613  }
614 
615  // for gray band
616  QDomElement channelElem = doc.createElement( QStringLiteral( "sld:GrayChannel" ) );
617  channelSelectionElem.appendChild( channelElem );
618 
619  // set band
620  QDomElement sourceChannelNameElem = doc.createElement( QStringLiteral( "sld:SourceChannelName" ) );
621  sourceChannelNameElem.appendChild( doc.createTextNode( QString::number( band() ) ) );
622  channelElem.appendChild( sourceChannelNameElem );
623  }
624 
625  // add ShadedRelief tag
626  QDomElement shadedReliefElem = doc.createElement( QStringLiteral( "sld:ShadedRelief" ) );
627  rasterSymbolizerElem.appendChild( shadedReliefElem );
628 
629  // brightnessOnly tag
630  QDomElement brightnessOnlyElem = doc.createElement( QStringLiteral( "sld:BrightnessOnly" ) );
631  brightnessOnlyElem.appendChild( doc.createTextNode( QStringLiteral( "true" ) ) );
632  shadedReliefElem.appendChild( brightnessOnlyElem );
633 
634  // ReliefFactor tag
635  QDomElement reliefFactorElem = doc.createElement( QStringLiteral( "sld:ReliefFactor" ) );
636  reliefFactorElem.appendChild( doc.createTextNode( QString::number( zFactor() ) ) );
637  shadedReliefElem.appendChild( reliefFactorElem );
638 
639  // altitude VendorOption tag
640  QDomElement altitudeVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
641  altitudeVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "altitude" ) );
642  altitudeVendorOptionElem.appendChild( doc.createTextNode( QString::number( altitude() ) ) );
643  shadedReliefElem.appendChild( altitudeVendorOptionElem );
644 
645  // azimuth VendorOption tag
646  QDomElement azimutVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
647  azimutVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "azimuth" ) );
648  azimutVendorOptionElem.appendChild( doc.createTextNode( QString::number( azimuth() ) ) );
649  shadedReliefElem.appendChild( azimutVendorOptionElem );
650 
651  // multidirectional VendorOption tag
652  QDomElement multidirectionalVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
653  multidirectionalVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "multidirectional" ) );
654  multidirectionalVendorOptionElem.appendChild( doc.createTextNode( QString::number( multiDirectional() ) ) );
655  shadedReliefElem.appendChild( multidirectionalVendorOptionElem );
656 }
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:58
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:265
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
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:587
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.
static const QRgb NODATA_COLOR
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:94
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:596
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:54
static QLatin1String LOGMESSAGE_TAG
OpenCL string for message logs.
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