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