QGIS API Documentation 3.37.0-Master (fa0c88a90b0)
qgselevationmap.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgselevationmap.cpp
3 --------------------------------------
4 Date : August 2022
5 Copyright : (C) 2022 by Martin Dobias
6 Email : wonder dot sk at gmail dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
16#include "qgselevationmap.h"
17
18#include "qgsrasterblock.h"
19
20#include <QPainter>
21#include <algorithm>
22#include <cmath>
23
24
25static const int ELEVATION_OFFSET = 7900;
26static const int ELEVATION_SCALE = 1000;
27
28
29QgsElevationMap::QgsElevationMap( const QSize &size, float devicePixelRatio )
30 : mElevationImage( size, QImage::Format_ARGB32 )
31{
32 mElevationImage.fill( 0 );
33 mElevationImage.setDevicePixelRatio( devicePixelRatio );
34}
35
37 : mElevationImage( image )
38{}
39
41 : mElevationImage( other.mElevationImage )
42{
43 mPainter.reset();
44}
45
47{
48 double zScaled = ( z + ELEVATION_OFFSET ) * ELEVATION_SCALE;
49 unsigned int zInt = static_cast<unsigned int>( std::clamp( zScaled, 0., 16777215. ) ); // make sure to fit into 3 bytes
50 return QRgb( zInt | ( static_cast< unsigned int >( 0xff ) << 24 ) );
51}
52
54{
55 unsigned int zScaled = colorRaw & 0xffffff;
56 return ( ( float ) zScaled ) / ELEVATION_SCALE - ELEVATION_OFFSET;
57}
58
59std::unique_ptr<QgsElevationMap> QgsElevationMap::fromRasterBlock( QgsRasterBlock *block )
60{
61 std::unique_ptr<QgsElevationMap> elevMap( new QgsElevationMap( QSize( block->width(), block->height() ) ) );
62 QRgb *dataPtr = reinterpret_cast<QRgb *>( elevMap->mElevationImage.bits() );
63 for ( int row = 0; row < block->height(); ++row )
64 {
65 for ( int col = 0; col < block->width(); ++col )
66 {
67 bool isNoData;
68 double value = block->valueAndNoData( row, col, isNoData );
69 if ( !isNoData )
70 *dataPtr = encodeElevation( static_cast<float>( value ) );
71 ++dataPtr;
72 }
73 }
74 return elevMap;
75}
76
78{
79 mPainter.reset();
80 mElevationImage = other.mElevationImage;
81 return *this;
82}
83
84void QgsElevationMap::applyEyeDomeLighting( QImage &img, int distance, float strength, float rendererScale ) const
85{
86 const int imgWidth = img.width(), imgHeight = img.height();
87 QRgb *imgPtr = reinterpret_cast<QRgb *>( img.bits() );
88 const QRgb *elevPtr = reinterpret_cast<const QRgb *>( mElevationImage.constBits() );
89
90 const int neighbours[] = { -1, 0, 1, 0, 0, -1, 0, 1 };
91 for ( int i = 0; i < imgWidth; ++i )
92 {
93 for ( int j = 0; j < imgHeight; ++j )
94 {
95 qgssize index = j * static_cast<qgssize>( imgWidth ) + i;
96 float factor = 0.0f;
97 float centerDepth = decodeElevation( elevPtr[ index ] );
98 if ( isNoData( elevPtr[ index ] ) )
99 continue;
100 for ( qgssize k = 0; k < 4; ++k )
101 {
102 float borderFactor = 1.0f;
103 int iNeighbour = std::clamp( i + distance * neighbours[2 * k], 0, imgWidth - 1 );
104 int jNeighbour = std::clamp( j + distance * neighbours[2 * k + 1], 0, imgHeight - 1 );
105 qgssize neighbourIndex = jNeighbour * static_cast<qgssize>( imgWidth ) + iNeighbour;
106
107 // neighbour is no data, we try to reach one pixel closer that is not no data
108 // and calculate a factor to take account of the reduction of the distance
109 if ( isNoData( elevPtr[ neighbourIndex ] ) )
110 {
111 if ( neighbours[2 * k] != 0 )
112 {
113 int corrI = iNeighbour;
114 int inc = neighbours[2 * k];
115 while ( corrI != i && isNoData( elevPtr[ neighbourIndex ] ) )
116 {
117 corrI -= inc;
118 neighbourIndex = jNeighbour * static_cast<qgssize>( imgWidth ) + corrI;
119 }
120 if ( corrI == i )
121 continue;
122
123 borderFactor = static_cast<float>( distance ) / static_cast<float>( std::abs( i - corrI ) ) ;
124 }
125
126 if ( neighbours[2 * k + 1] != 0 )
127 {
128 int corrJ = jNeighbour;
129 int inc = neighbours[2 * k + 1];
130 while ( corrJ != j && isNoData( elevPtr[ neighbourIndex ] ) )
131 {
132 corrJ -= inc;
133 neighbourIndex = corrJ * static_cast<qgssize>( imgWidth ) + iNeighbour;
134 }
135 if ( corrJ == j )
136 continue;
137
138 borderFactor = static_cast<float>( distance ) / static_cast<float>( std::abs( j - corrJ ) ) ;
139 }
140 }
141
142 // neighbour is outside the extent (right or left), we take the pixel on border
143 // and calculate a factor to take account of the reduction of the distance
144 if ( neighbours[2 * k] != 0 && std::abs( iNeighbour - i ) < distance )
145 {
146 if ( i > iNeighbour )
147 borderFactor = static_cast<float>( distance ) / static_cast<float>( i ) ;
148 else if ( i < iNeighbour )
149 borderFactor = static_cast<float>( distance ) / static_cast<float>( ( imgWidth - i - 1 ) ) ;
150 }
151 if ( neighbours[2 * k + 1] != 0 && std::abs( jNeighbour - j ) < distance )
152 {
153 if ( j > jNeighbour )
154 borderFactor = static_cast<float>( distance ) / static_cast<float>( j ) ;
155 else if ( j < jNeighbour )
156 borderFactor = static_cast<float>( distance ) / static_cast<float>( imgHeight - j - 1 ) ;
157 }
158
159 float neighbourDepth = decodeElevation( elevPtr[ neighbourIndex ] );
160 factor += std::max<float>( 0, -( centerDepth - neighbourDepth ) * borderFactor );
161 }
162 float shade = expf( -factor * strength / rendererScale );
163
164 QRgb c = imgPtr[ index ];
165 imgPtr[ index ] = qRgba( static_cast<int>( static_cast<float>( qRed( c ) ) * shade ),
166 static_cast<int>( static_cast<float>( qGreen( c ) ) * shade ),
167 static_cast<int>( static_cast<float>( qBlue( c ) ) * shade )
168 , qAlpha( c ) );
169 }
170 }
171}
172
173void QgsElevationMap::applyHillshading( QImage &img, bool multiDirectional, double altitude, double azimuth, double zFactor, double cellSizeX, double cellSizeY ) const
174{
175 // algs from src/raster/qgshillshaderenderer.cpp
176 double altRad = altitude * M_PI / 180.0;
177 double cos_altRadian = std::cos( altRad );
178 double sin_altRadian = std::sin( altRad );
179 double cos_alt_mul_z = cos_altRadian * zFactor ;
180 double azimuthRad = -1 * azimuth * M_PI / 180.0;
181 double cos_az_mul_cos_alt_mul_z = std::cos( azimuthRad ) * cos_alt_mul_z;
182 double sin_az_mul_cos_alt_mul_z = std::sin( azimuthRad ) * cos_alt_mul_z;
183
184 // From qgshillshaderenderer.cpp, -32.87001872802012 comes from ( 127.0 * std::cos(225.0 * M_PI / 180.0))
185 // but this expression actually equals -89.8025612106916. To be consistent, we keep -32.87001872802012
186 // and divide it by 254 because, here, values are between 0.0 and 1.0, not between 0 and 255
187 double cos225_az_mul_cos_alt_mul_z_mul_0_5 = -32.87001872802012 * cos_alt_mul_z / 255;
188 double square_z = zFactor * zFactor;
189
190 QRgb *imgPtr = reinterpret_cast<QRgb *>( img.bits() );
191 const QRgb *elevPtr = reinterpret_cast<const QRgb *>( mElevationImage.constBits() );
192
193 const int imgWidth = img.width(), imgHeight = img.height();
194
195 auto colRowToIndex = [&]( int c, int r )->qgssize
196 {
197 return static_cast<qgssize>( r ) * imgWidth + c ;
198 };
199
200 float noData = noDataValue();
201
202 // Elevation value matrix
203 // 0 1 2 11 12 13
204 // 3 4 5 21 22 23
205 // 6 7 8 31 32 33
206 float pixelValues[9];
207
208 for ( int rowC = 0; rowC < imgHeight ; ++rowC )
209 {
210 int rowU = std::max( 0, rowC - 1 );
211 int rowD = std::min( imgHeight - 1, rowC + 1 );
212
213 pixelValues[1] = decodeElevation( elevPtr[colRowToIndex( 0, rowU )] );
214 pixelValues[4] = decodeElevation( elevPtr[colRowToIndex( 0, rowC )] );
215 pixelValues[7] = decodeElevation( elevPtr[colRowToIndex( 0, rowD )] );
216
217 pixelValues[2] = decodeElevation( elevPtr[colRowToIndex( 1, rowU )] );
218 pixelValues[5] = decodeElevation( elevPtr[colRowToIndex( 1, rowC )] );
219 pixelValues[8] = decodeElevation( elevPtr[colRowToIndex( 1, rowD )] );
220
221
222 for ( int colC = 0; colC < imgWidth ; ++colC )
223 {
224 qgssize centerIndex = colRowToIndex( colC, rowC );
225 int colR = std::min( imgWidth - 1, colC + 1 );
226
227 pixelValues[0] = pixelValues[1];
228 pixelValues[3] = pixelValues[4];
229 pixelValues[6] = pixelValues[7];
230
231 pixelValues[1] = pixelValues[2];
232 pixelValues[4] = pixelValues[5];
233 pixelValues[7] = pixelValues[8];
234
235 pixelValues[2] = decodeElevation( elevPtr[colRowToIndex( colR, rowU )] );
236 pixelValues[5] = decodeElevation( elevPtr[colRowToIndex( colR, rowC )] );
237 pixelValues[8] = decodeElevation( elevPtr[colRowToIndex( colR, rowD )] );
238
239 if ( elevPtr[centerIndex] != 0 )
240 {
241 // This is center cell. Use this in place of nodata neighbors
242 const float x22 = pixelValues[4];
243 if ( x22 == noData )
244 continue;
245
246 const float x11 = ( pixelValues[0] == noData ) ? x22 : pixelValues[0];
247 const float x21 = ( pixelValues[3] == noData ) ? x22 : pixelValues[3];
248 const float x31 = ( pixelValues[6] == noData ) ? x22 : pixelValues[6];
249
250 const float x12 = ( pixelValues[1] == noData ) ? x22 : pixelValues[1];
251 const float x32 = ( pixelValues[7] == noData ) ? x22 : pixelValues[7];
252
253 const float x13 = ( pixelValues[2] == noData ) ? x22 : pixelValues[2];
254 const float x23 = ( pixelValues[5] == noData ) ? x22 : pixelValues[5];
255 const float x33 = ( pixelValues[8] == noData ) ? x22 : pixelValues[8];
256
257 const double derX = static_cast<double>( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellSizeX );
258 const double derY = static_cast<double>( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellSizeY );
259
260 double shade = 0;
261 if ( !multiDirectional )
262 {
263 // Standard single direction hillshade
264 shade = std::clamp( ( sin_altRadian -
265 ( derY * cos_az_mul_cos_alt_mul_z -
266 derX * sin_az_mul_cos_alt_mul_z ) ) /
267 std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) ),
268 0.0, 1.0 );
269 }
270 else
271 {
272 // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
273 // Fast formula from GDAL DEM
274 const double xx = derX * derX;
275 const double yy = derY * derY;
276 const double xx_plus_yy = xx + yy;
277 // Flat?
278 if ( xx_plus_yy == 0.0 )
279 {
280 shade = std::clamp( sin_altRadian, 0.0, 1.0 );
281 }
282 else
283 {
284 // ... then the shade value from different azimuth
285 double val225_mul_0_5 = sin_altRadian * 0.5 +
286 ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_0_5;
287 val225_mul_0_5 = ( val225_mul_0_5 <= 0.0 ) ? 0.0 : val225_mul_0_5;
288 double val270_mul_0_5 = sin_altRadian * 0.5 -
289 derX * cos_alt_mul_z * 0.5;
290 val270_mul_0_5 = ( val270_mul_0_5 <= 0.0 ) ? 0.0 : val270_mul_0_5;
291 double val315_mul_0_5 = sin_altRadian * 0.5 +
292 ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_0_5;
293 val315_mul_0_5 = ( val315_mul_0_5 <= 0.0 ) ? 0.0 : val315_mul_0_5;
294 double val360_mul_0_5 = sin_altRadian * 0.5 -
295 derY * cos_alt_mul_z * 0.5;
296 val360_mul_0_5 = ( val360_mul_0_5 <= 0.0 ) ? 0.0 : val360_mul_0_5;
297
298 // ... then the weighted shading
299 const double weight_225 = 0.5 * xx_plus_yy - derX * derY;
300 const double weight_270 = xx;
301 const double weight_315 = xx_plus_yy - weight_225;
302 const double weight_360 = yy;
303 const double cang_mul_127 = (
304 ( weight_225 * val225_mul_0_5 +
305 weight_270 * val270_mul_0_5 +
306 weight_315 * val315_mul_0_5 +
307 weight_360 * val360_mul_0_5 ) / xx_plus_yy ) /
308 ( 1 + square_z * xx_plus_yy );
309
310 shade = std::clamp( cang_mul_127, 0.0, 1.0 );
311 }
312 }
313
314 QRgb c = imgPtr[ centerIndex ];
315 imgPtr[ centerIndex ] = qRgba( static_cast<int>( static_cast<float>( qRed( c ) ) * shade ),
316 static_cast<int>( static_cast<float>( qGreen( c ) ) * shade ),
317 static_cast<int>( static_cast<float>( qBlue( c ) ) * shade )
318 , qAlpha( c ) );
319 }
320 }
321 }
322}
323
325{
326 if ( !mPainter )
327 {
328 mPainter.reset( new QPainter );
329 mPainter->begin( &mElevationImage );
330 }
331 return mPainter.get();
332
333}
334
336{
337 if ( otherElevationMap.mElevationImage.size() != mElevationImage.size() )
338 {
339 QgsDebugMsgLevel( QStringLiteral( "Elevation map with different sizes can not be combined" ), 4 );
340 return;
341 }
342
343 QRgb *elevPtr = reinterpret_cast<QRgb *>( mElevationImage.bits() );
344 const QRgb *otherElevPtr = reinterpret_cast<const QRgb *>( otherElevationMap.mElevationImage.constBits() );
345
346 int width = mElevationImage.width();
347 int height = mElevationImage.height();
348
349 for ( int row = 0; row < height; ++row )
350 {
351 for ( int col = 0; col < width; ++col )
352 {
353 qgssize index = col + static_cast<qgssize>( row ) * width ;
354 if ( !isNoData( otherElevPtr[index] ) )
355 {
356 switch ( method )
357 {
359 {
360 double elev = decodeElevation( elevPtr[index] );
361 double otherElev = decodeElevation( otherElevPtr[index] );
362 if ( otherElev > elev )
363 elevPtr[index] = otherElevPtr[index];
364 }
365 break;
367 elevPtr[index] = otherElevPtr[index];
368 break;
369 }
370 }
371 }
372 }
373}
374
375
377{
378 return !mElevationImage.isNull();
379}
380
381void QgsElevationMap::fillWithRasterBlock( QgsRasterBlock *block, int top, int left, double zScale, double offset )
382{
383 QRgb *dataPtr = reinterpret_cast<QRgb *>( mElevationImage.bits() );
384
385 int widthMap = mElevationImage.width();
386 int heightMap = mElevationImage.height();
387 int widthBlock = block->width();
388 int combinedHeight = std::min( mElevationImage.height(), block->height() + top );
389 int combinedWidth = std::min( widthMap, widthBlock + left );
390 for ( int row = std::max( 0, top ); row < combinedHeight; ++row )
391 {
392 if ( row >= heightMap )
393 continue;
394
395 for ( int col = std::max( 0, left ); col < combinedWidth; ++col )
396 {
397 if ( col >= widthMap )
398 continue;
399
400 bool isNoData = true;
401 double value = block->valueAndNoData( ( row - top ), ( col - left ), isNoData );
402 qgssize index = col + static_cast<qgssize>( row ) * widthMap;
403 if ( isNoData )
404 dataPtr[index] = 0;
405 else
406 dataPtr[index] = encodeElevation( static_cast<float>( value * zScale + offset ) );
407 }
408 }
409}
ElevationMapCombineMethod
Methods used to select the elevation when two elevation maps are combined.
Definition: qgis.h:3966
@ NewerElevation
Keep the new elevation regardless of its value if it is not null.
@ HighestElevation
Keep the highest elevation if it is not null.
Stores digital elevation model in a raster image which may get updated as a part of map layer renderi...
void combine(const QgsElevationMap &otherElevationMap, Qgis::ElevationMapCombineMethod method)
Combines this elevation map with otherElevationMap.
static QRgb encodeElevation(float z)
Converts elevation value to an actual color.
QgsElevationMap()=default
Default constructor.
void applyEyeDomeLighting(QImage &image, int distance, float strength, float rendererScale) const
Applies eye dome lighting effect to the given image.
static std::unique_ptr< QgsElevationMap > fromRasterBlock(QgsRasterBlock *block)
Creates an elevation map based on data from the given raster block.
QPainter * painter() const
Returns painter to the underlying QImage with elevations.
float noDataValue() const
Returns the no data value for the elevation map.
bool isNoData(QRgb colorRaw) const
Returns whether the encoded value is a no data value.
bool isValid() const
Returns whether the elevation map is valid.
static float decodeElevation(QRgb colorRaw)
Converts a color back to elevation value.
void fillWithRasterBlock(QgsRasterBlock *block, int top, int left, double zScale=1.0, double offset=0.0)
Fills the elevation map with values contains in a raster block starting from position defined by top ...
void applyHillshading(QImage &image, bool multiDirectional, double altitude, double azimuth, double zFactor, double cellSizeX, double cellSizeY) const
Applies hill shading effect to the given image.
QgsElevationMap & operator=(const QgsElevationMap &other)
Raster data container.
int height() const
Returns the height (number of rows) of the raster block.
double valueAndNoData(int row, int column, bool &isNoData) const
Reads a single value from the pixel at row and column, if type of block is numeric.
int width() const
Returns the width (number of columns) of the raster block.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
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:5686
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39