Branch data Line data Source code
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 : 0 : QgsHillshadeRenderer::QgsHillshadeRenderer( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ):
38 : 0 : QgsRasterRenderer( input, QStringLiteral( "hillshade" ) )
39 : 0 : , mBand( band )
40 : 0 : , mZFactor( 1 )
41 : 0 : , mLightAngle( lightAngle )
42 : 0 : , mLightAzimuth( lightAzimuth )
43 : 0 : , mMultiDirectional( false )
44 : 0 : {
45 : :
46 : 0 : }
47 : :
48 : 0 : QgsHillshadeRenderer *QgsHillshadeRenderer::clone() const
49 : : {
50 : 0 : QgsHillshadeRenderer *r = new QgsHillshadeRenderer( nullptr, mBand, mLightAzimuth, mLightAngle );
51 : 0 : r->copyCommonProperties( this );
52 : :
53 : 0 : r->setZFactor( mZFactor );
54 : 0 : r->setMultiDirectional( mMultiDirectional );
55 : 0 : return r;
56 : 0 : }
57 : :
58 : 0 : QgsRasterRenderer *QgsHillshadeRenderer::create( const QDomElement &elem, QgsRasterInterface *input )
59 : : {
60 : 0 : if ( elem.isNull() )
61 : : {
62 : 0 : return nullptr;
63 : : }
64 : :
65 : 0 : int band = elem.attribute( QStringLiteral( "band" ), QStringLiteral( "0" ) ).toInt();
66 : 0 : double azimuth = elem.attribute( QStringLiteral( "azimuth" ), QStringLiteral( "315" ) ).toDouble();
67 : 0 : double angle = elem.attribute( QStringLiteral( "angle" ), QStringLiteral( "45" ) ).toDouble();
68 : 0 : double zFactor = elem.attribute( QStringLiteral( "zfactor" ), QStringLiteral( "1" ) ).toDouble();
69 : 0 : bool multiDirectional = elem.attribute( QStringLiteral( "multidirection" ), QStringLiteral( "0" ) ).toInt();
70 : 0 : QgsHillshadeRenderer *r = new QgsHillshadeRenderer( input, band, azimuth, angle );
71 : 0 : r->readXml( elem );
72 : :
73 : 0 : r->setZFactor( zFactor );
74 : 0 : r->setMultiDirectional( multiDirectional );
75 : 0 : return r;
76 : 0 : }
77 : :
78 : 0 : void QgsHillshadeRenderer::writeXml( QDomDocument &doc, QDomElement &parentElem ) const
79 : : {
80 : 0 : if ( parentElem.isNull() )
81 : : {
82 : 0 : return;
83 : : }
84 : :
85 : 0 : QDomElement rasterRendererElem = doc.createElement( QStringLiteral( "rasterrenderer" ) );
86 : 0 : _writeXml( doc, rasterRendererElem );
87 : :
88 : 0 : rasterRendererElem.setAttribute( QStringLiteral( "band" ), mBand );
89 : 0 : rasterRendererElem.setAttribute( QStringLiteral( "azimuth" ), QString::number( mLightAzimuth ) );
90 : 0 : rasterRendererElem.setAttribute( QStringLiteral( "angle" ), QString::number( mLightAngle ) );
91 : 0 : rasterRendererElem.setAttribute( QStringLiteral( "zfactor" ), QString::number( mZFactor ) );
92 : 0 : rasterRendererElem.setAttribute( QStringLiteral( "multidirection" ), QString::number( mMultiDirectional ) );
93 : 0 : parentElem.appendChild( rasterRendererElem );
94 : 0 : }
95 : :
96 : 0 : QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
97 : : {
98 : : Q_UNUSED( bandNo )
99 : 0 : std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock() );
100 : 0 : if ( !mInput )
101 : : {
102 : 0 : QgsDebugMsg( QStringLiteral( "No input raster!" ) );
103 : 0 : return outputBlock.release();
104 : : }
105 : :
106 : 0 : std::shared_ptr< QgsRasterBlock > inputBlock( mInput->block( mBand, extent, width, height, feedback ) );
107 : :
108 : 0 : if ( !inputBlock || inputBlock->isEmpty() )
109 : : {
110 : 0 : QgsDebugMsg( QStringLiteral( "No raster data!" ) );
111 : 0 : return outputBlock.release();
112 : : }
113 : :
114 : 0 : std::shared_ptr< QgsRasterBlock > alphaBlock;
115 : :
116 : 0 : if ( mAlphaBand > 0 && mBand != mAlphaBand )
117 : : {
118 : 0 : alphaBlock.reset( mInput->block( mAlphaBand, extent, width, height, feedback ) );
119 : 0 : if ( !alphaBlock || alphaBlock->isEmpty() )
120 : : {
121 : : // TODO: better to render without alpha
122 : 0 : return outputBlock.release();
123 : : }
124 : 0 : }
125 : 0 : else if ( mAlphaBand > 0 )
126 : : {
127 : 0 : alphaBlock = inputBlock;
128 : 0 : }
129 : :
130 : 0 : if ( !outputBlock->reset( Qgis::ARGB32_Premultiplied, width, height ) )
131 : : {
132 : 0 : return outputBlock.release();
133 : : }
134 : :
135 : : // Starting the computation
136 : :
137 : : // Common pre-calculated values
138 : 0 : float cellXSize = static_cast<float>( extent.width() ) / width;
139 : 0 : float cellYSize = static_cast<float>( extent.height() ) / height;
140 : 0 : float zenithRad = static_cast<float>( std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0 );
141 : 0 : float azimuthRad = static_cast<float>( -1 * mLightAzimuth * M_PI / 180.0 );
142 : 0 : float cosZenithRad = std::cos( zenithRad );
143 : 0 : float sinZenithRad = std::sin( zenithRad );
144 : :
145 : : // For fast formula from GDAL DEM
146 : 0 : float cos_alt_mul_z = cosZenithRad * static_cast<float>( mZFactor );
147 : 0 : float cos_az_mul_cos_alt_mul_z = std::cos( azimuthRad ) * cos_alt_mul_z;
148 : 0 : float sin_az_mul_cos_alt_mul_z = std::sin( azimuthRad ) * cos_alt_mul_z;
149 : 0 : float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0f * cos_az_mul_cos_alt_mul_z;
150 : 0 : float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0f * sin_az_mul_cos_alt_mul_z;
151 : 0 : float square_z = static_cast<float>( mZFactor * mZFactor );
152 : 0 : float sin_altRadians_mul_254 = 254.0f * sinZenithRad;
153 : :
154 : : // For multi directional
155 : 0 : float sin_altRadians_mul_127 = 127.0f * sinZenithRad;
156 : : // 127.0 * std::cos(225.0 * M_PI / 180.0) = -32.87001872802012
157 : 0 : float cos225_az_mul_cos_alt_mul_z_mul_127 = -32.87001872802012f * cos_alt_mul_z;
158 : 0 : float cos_alt_mul_z_mul_127 = 127.0f * cos_alt_mul_z;
159 : :
160 : 0 : 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()
166 : : && QgsOpenClUtils::available()
167 : : && ( ! mRasterTransparency || mRasterTransparency->isEmpty() )
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( QLatin1String( "__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 = std::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( ) ) ),
360 : : QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
361 : : QgsOpenClUtils::setEnabled( false );
362 : : QgsMessageLog::logMessage( QObject::tr( "OpenCL has been disabled, you can re-enable it in the options dialog." ),
363 : : QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
364 : : }
365 : :
366 : : } // End of OpenCL processing path
367 : : else // Use the CPU and the original algorithm
368 : : {
369 : :
370 : : #endif
371 : :
372 : 0 : for ( qgssize i = 0; i < static_cast<qgssize>( height ); i++ )
373 : : {
374 : :
375 : 0 : for ( qgssize j = 0; j < static_cast<qgssize>( width ); j++ )
376 : : {
377 : :
378 : 0 : if ( inputBlock->isNoData( i, j ) )
379 : : {
380 : 0 : outputBlock->setColor( static_cast<int>( i ), static_cast<int>( j ), defaultNodataColor );
381 : 0 : continue;
382 : : }
383 : :
384 : : qgssize iUp, iDown, jLeft, jRight;
385 : 0 : if ( i == 0 )
386 : : {
387 : 0 : iUp = i;
388 : 0 : iDown = i + 1;
389 : 0 : }
390 : 0 : else if ( i < static_cast<qgssize>( height ) - 1 )
391 : : {
392 : 0 : iUp = i - 1;
393 : 0 : iDown = i + 1;
394 : 0 : }
395 : : else
396 : : {
397 : 0 : iUp = i - 1;
398 : 0 : iDown = i;
399 : : }
400 : :
401 : 0 : if ( j == 0 )
402 : : {
403 : 0 : jLeft = j;
404 : 0 : jRight = j + 1;
405 : 0 : }
406 : 0 : else if ( j < static_cast<qgssize>( width ) - 1 )
407 : : {
408 : 0 : jLeft = j - 1;
409 : 0 : jRight = j + 1;
410 : 0 : }
411 : : else
412 : : {
413 : 0 : jLeft = j - 1;
414 : 0 : 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 : 0 : x22 = inputBlock->value( i, j );
429 : :
430 : 0 : x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft );
431 : 0 : x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft );
432 : 0 : x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft );
433 : :
434 : 0 : x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j );
435 : : // x22
436 : 0 : x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j );
437 : :
438 : 0 : x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight );
439 : 0 : x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight );
440 : 0 : x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight );
441 : :
442 : 0 : double derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize );
443 : 0 : double derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize );
444 : :
445 : : // Fast formula
446 : :
447 : : double grayValue;
448 : 0 : if ( !mMultiDirectional )
449 : : {
450 : : // Standard single direction hillshade
451 : 0 : grayValue = std::clamp( ( sin_altRadians_mul_254 -
452 : 0 : ( derY * cos_az_mul_cos_alt_mul_z_mul_254 -
453 : 0 : derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) /
454 : 0 : std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) ),
455 : 0 : 0.0, 255.0 );
456 : 0 : }
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 : 0 : const float xx = derX * derX;
462 : 0 : const float yy = derY * derY;
463 : 0 : const float xx_plus_yy = xx + yy;
464 : : // Flat?
465 : 0 : if ( xx_plus_yy == 0.0 )
466 : : {
467 : 0 : grayValue = std::clamp( static_cast<float>( 1.0 + sin_altRadians_mul_254 ), 0.0f, 255.0f );
468 : 0 : }
469 : : else
470 : : {
471 : : // ... then the shade value from different azimuth
472 : 0 : float val225_mul_127 = sin_altRadians_mul_127 +
473 : 0 : ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
474 : 0 : val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
475 : 0 : float val270_mul_127 = sin_altRadians_mul_127 -
476 : 0 : derX * cos_alt_mul_z_mul_127;
477 : 0 : val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
478 : 0 : float val315_mul_127 = sin_altRadians_mul_127 +
479 : 0 : ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
480 : 0 : val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
481 : 0 : float val360_mul_127 = sin_altRadians_mul_127 -
482 : 0 : derY * cos_alt_mul_z_mul_127;
483 : 0 : val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
484 : :
485 : : // ... then the weighted shading
486 : 0 : const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
487 : 0 : const float weight_270 = xx;
488 : 0 : const float weight_315 = xx_plus_yy - weight_225;
489 : 0 : const float weight_360 = yy;
490 : 0 : const float cang_mul_127 = (
491 : 0 : ( weight_225 * val225_mul_127 +
492 : 0 : weight_270 * val270_mul_127 +
493 : 0 : weight_315 * val315_mul_127 +
494 : 0 : weight_360 * val360_mul_127 ) / xx_plus_yy ) /
495 : 0 : ( 1 + square_z * xx_plus_yy );
496 : :
497 : 0 : grayValue = std::clamp( 1.0f + cang_mul_127, 0.0f, 255.0f );
498 : : }
499 : : }
500 : :
501 : 0 : double currentAlpha = mOpacity;
502 : 0 : if ( mRasterTransparency )
503 : : {
504 : 0 : currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0;
505 : 0 : }
506 : 0 : if ( mAlphaBand > 0 )
507 : : {
508 : 0 : currentAlpha *= alphaBlock->value( i ) / 255.0;
509 : 0 : }
510 : :
511 : 0 : if ( qgsDoubleNear( currentAlpha, 1.0 ) )
512 : : {
513 : 0 : outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) );
514 : 0 : }
515 : : else
516 : : {
517 : 0 : outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
518 : : }
519 : 0 : }
520 : 0 : }
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 : 0 : return outputBlock.release();
539 : 0 : }
540 : :
541 : 0 : QList<int> QgsHillshadeRenderer::usesBands() const
542 : : {
543 : 0 : QList<int> bandList;
544 : 0 : if ( mBand != -1 )
545 : : {
546 : 0 : bandList << mBand;
547 : 0 : }
548 : 0 : return bandList;
549 : :
550 : 0 : }
551 : :
552 : 0 : void QgsHillshadeRenderer::setBand( int bandNo )
553 : : {
554 : 0 : if ( bandNo > mInput->bandCount() || bandNo <= 0 )
555 : : {
556 : 0 : return;
557 : : }
558 : 0 : mBand = bandNo;
559 : 0 : }
560 : :
561 : 0 : 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 : 0 : return ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellsize );
567 : : }
568 : :
569 : 0 : 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 : 0 : return ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellsize );
575 : : }
576 : :
577 : 0 : void QgsHillshadeRenderer::toSld( QDomDocument &doc, QDomElement &element, const QVariantMap &props ) const
578 : : {
579 : : // create base structure
580 : 0 : QgsRasterRenderer::toSld( doc, element, props );
581 : :
582 : : // look for RasterSymbolizer tag
583 : 0 : QDomNodeList elements = element.elementsByTagName( QStringLiteral( "sld:RasterSymbolizer" ) );
584 : 0 : if ( elements.size() == 0 )
585 : 0 : return;
586 : :
587 : : // there SHOULD be only one
588 : 0 : 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 : 0 : if ( band() != 1 )
594 : : {
595 : 0 : QDomElement channelSelectionElem = doc.createElement( QStringLiteral( "sld:ChannelSelection" ) );
596 : 0 : elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral( "sld:Opacity" ) );
597 : 0 : if ( elements.size() != 0 )
598 : : {
599 : 0 : rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
600 : 0 : }
601 : : else
602 : : {
603 : 0 : elements = rasterSymbolizerElem.elementsByTagName( QStringLiteral( "sld:Geometry" ) );
604 : 0 : if ( elements.size() != 0 )
605 : : {
606 : 0 : rasterSymbolizerElem.insertAfter( channelSelectionElem, elements.at( 0 ) );
607 : 0 : }
608 : : else
609 : : {
610 : 0 : rasterSymbolizerElem.insertBefore( channelSelectionElem, rasterSymbolizerElem.firstChild() );
611 : : }
612 : : }
613 : :
614 : : // for gray band
615 : 0 : QDomElement channelElem = doc.createElement( QStringLiteral( "sld:GrayChannel" ) );
616 : 0 : channelSelectionElem.appendChild( channelElem );
617 : :
618 : : // set band
619 : 0 : QDomElement sourceChannelNameElem = doc.createElement( QStringLiteral( "sld:SourceChannelName" ) );
620 : 0 : sourceChannelNameElem.appendChild( doc.createTextNode( QString::number( band() ) ) );
621 : 0 : channelElem.appendChild( sourceChannelNameElem );
622 : 0 : }
623 : :
624 : : // add ShadedRelief tag
625 : 0 : QDomElement shadedReliefElem = doc.createElement( QStringLiteral( "sld:ShadedRelief" ) );
626 : 0 : rasterSymbolizerElem.appendChild( shadedReliefElem );
627 : :
628 : : // brightnessOnly tag
629 : 0 : QDomElement brightnessOnlyElem = doc.createElement( QStringLiteral( "sld:BrightnessOnly" ) );
630 : 0 : brightnessOnlyElem.appendChild( doc.createTextNode( QStringLiteral( "true" ) ) );
631 : 0 : shadedReliefElem.appendChild( brightnessOnlyElem );
632 : :
633 : : // ReliefFactor tag
634 : 0 : QDomElement reliefFactorElem = doc.createElement( QStringLiteral( "sld:ReliefFactor" ) );
635 : 0 : reliefFactorElem.appendChild( doc.createTextNode( QString::number( zFactor() ) ) );
636 : 0 : shadedReliefElem.appendChild( reliefFactorElem );
637 : :
638 : : // altitude VendorOption tag
639 : 0 : QDomElement altitudeVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
640 : 0 : altitudeVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "altitude" ) );
641 : 0 : altitudeVendorOptionElem.appendChild( doc.createTextNode( QString::number( altitude() ) ) );
642 : 0 : shadedReliefElem.appendChild( altitudeVendorOptionElem );
643 : :
644 : : // azimuth VendorOption tag
645 : 0 : QDomElement azimutVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
646 : 0 : azimutVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "azimuth" ) );
647 : 0 : azimutVendorOptionElem.appendChild( doc.createTextNode( QString::number( azimuth() ) ) );
648 : 0 : shadedReliefElem.appendChild( azimutVendorOptionElem );
649 : :
650 : : // multidirectional VendorOption tag
651 : 0 : QDomElement multidirectionalVendorOptionElem = doc.createElement( QStringLiteral( "sld:VendorOption" ) );
652 : 0 : multidirectionalVendorOptionElem.setAttribute( QStringLiteral( "name" ), QStringLiteral( "multidirectional" ) );
653 : 0 : multidirectionalVendorOptionElem.appendChild( doc.createTextNode( QString::number( multiDirectional() ) ) );
654 : 0 : shadedReliefElem.appendChild( multidirectionalVendorOptionElem );
655 : 0 : }
|
