22
22
#include " qgsrasterinterface.h"
23
23
#include " qgsrasterblock.h"
24
24
#include " qgsrectangle.h"
25
+ #include " qgssettings.h"
25
26
#include < memory>
26
27
28
+ #ifdef QGISDEBUG
29
+ #include " qgsmessagelog.h"
30
+ #include < chrono>
31
+ #endif
32
+
27
33
QgsHillshadeRenderer::QgsHillshadeRenderer ( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ):
28
34
QgsRasterRenderer( input, QStringLiteral( " hillshade" ) )
29
35
, mBand( band )
@@ -122,21 +128,34 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
122
128
return outputBlock.release ();
123
129
}
124
130
125
- double cellXSize = extent.width () / double ( width );
126
- double cellYSize = extent.height () / double ( height );
127
- double zenithRad = std::max ( 0.0 , 90 - mLightAngle ) * M_PI / 180.0 ;
128
- double azimuthRad = -1 * mLightAzimuth * M_PI / 180.0 ;
129
- double cosZenithRad = std::cos ( zenithRad );
130
- double sinZenithRad = std::sin ( zenithRad );
131
-
132
- // Multi direction hillshade: http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
133
- double angle0Rad = ( -1 * mLightAzimuth - 45 - 45 * 0.5 ) * M_PI / 180.0 ;
134
- double angle1Rad = ( -1 * mLightAzimuth - 45 * 0.5 ) * M_PI / 180.0 ;
135
- double angle2Rad = ( -1 * mLightAzimuth + 45 * 0.5 ) * M_PI / 180.0 ;
136
- double angle3Rad = ( -1 * mLightAzimuth + 45 + 45 * 0.5 ) * M_PI / 180.0 ;
131
+ float cellXSize = extent.width () / static_cast <float >( width );
132
+ float cellYSize = extent.height () / static_cast <float >( height );
133
+ float zenithRad = std::max ( 0.0 , 90 - mLightAngle ) * M_PI / 180.0 ;
134
+ float azimuthRad = -1 * mLightAzimuth * M_PI / 180.0 ;
135
+ float cosZenithRad = std::cos ( zenithRad );
136
+ float sinZenithRad = std::sin ( zenithRad );
137
+
138
+ // For fast formula from GDAL DEM
139
+ float cos_alt_mul_z = cosZenithRad * mZFactor ;
140
+ float cos_az_mul_cos_alt_mul_z = std::cos ( azimuthRad ) * cos_alt_mul_z;
141
+ float sin_az_mul_cos_alt_mul_z = std::sin ( azimuthRad ) * cos_alt_mul_z;
142
+ float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0 * cos_az_mul_cos_alt_mul_z;
143
+ float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0 * sin_az_mul_cos_alt_mul_z;
144
+ float square_z = mZFactor * mZFactor ;
145
+ float sinZenithRad_mul_254 = 254.0 * sinZenithRad;
146
+
147
+ // For multi directional
148
+ float sinZenithRad_mul_127 = 127.0 * sinZenithRad;
149
+ // 127.0 * std::cos(225.0 * M_PI / 180.0) = -32.87001872802012
150
+ float cos225_az_mul_cos_alt_mul_z_mul_127 = -32 .87001872802012f * cos_alt_mul_z;
151
+ float cos_alt_mul_z_mul_127 = 127.0 * cos_alt_mul_z;
137
152
138
153
QRgb myDefaultColor = NODATA_COLOR;
139
154
155
+ #ifdef QGISDEBUG
156
+ std::chrono::time_point<std::chrono::system_clock> startTime ( std::chrono::system_clock::now () );
157
+ #endif
158
+
140
159
for ( qgssize i = 0 ; i < ( qgssize )height; i++ )
141
160
{
142
161
@@ -182,15 +201,15 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
182
201
jRight = j;
183
202
}
184
203
185
- double x11;
186
- double x21;
187
- double x31;
188
- double x12;
189
- double x22; // Working cell
190
- double x32;
191
- double x13;
192
- double x23;
193
- double x33;
204
+ float x11;
205
+ float x21;
206
+ float x31;
207
+ float x12;
208
+ float x22; // Working cell
209
+ float x32;
210
+ float x13;
211
+ float x23;
212
+ float x33;
194
213
195
214
// This is center cell. It is not nodata. Use this in place of nodata neighbors
196
215
x22 = inputBlock->value ( i, j );
@@ -207,43 +226,65 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
207
226
x23 = inputBlock->isNoData ( i, jRight ) ? x22 : inputBlock->value ( i, jRight );
208
227
x33 = inputBlock->isNoData ( iDown, jRight ) ? x22 : inputBlock->value ( iDown, jRight );
209
228
210
- double derX = calcFirstDerX ( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize );
211
- double derY = calcFirstDerY ( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize );
212
-
213
- double slopeRad = std::atan ( mZFactor * std::sqrt ( derX * derX + derY * derY ) );
214
- double aspectRad = std::atan2 ( derX, -derY );
229
+ float derX = calcFirstDerX ( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize );
230
+ float derY = calcFirstDerY ( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize );
215
231
216
-
217
- double grayValue;
232
+ float grayValue;
218
233
if ( !mMultiDirectional )
219
234
{
220
235
// Standard single direction hillshade
221
- grayValue = qBound ( 0.0 , 255.0 * ( cosZenithRad * std::cos ( slopeRad )
222
- + sinZenithRad * std::sin ( slopeRad )
223
- * std::cos ( azimuthRad - aspectRad ) ), 255.0 );
236
+ // Fast formula from GDAL DEM
237
+ grayValue = qBound ( 0 .0f , ( sinZenithRad_mul_254 -
238
+ ( derY * cos_az_mul_cos_alt_mul_z_mul_254 -
239
+ derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) /
240
+ std::sqrt ( 1 + square_z * ( derX * derX + derY * derY ) )
241
+ , 255 .0f );
224
242
}
225
243
else
226
244
{
227
245
// Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
228
- double weight0 = std::sin ( aspectRad - angle0Rad );
229
- double weight1 = std::sin ( aspectRad - angle1Rad );
230
- double weight2 = std::sin ( aspectRad - angle2Rad );
231
- double weight3 = std::sin ( aspectRad - angle3Rad );
232
- weight0 = weight0 * weight0;
233
- weight1 = weight1 * weight1;
234
- weight2 = weight2 * weight2;
235
- weight3 = weight3 * weight3;
236
-
237
- double cosSlope = cosZenithRad * std::cos ( slopeRad );
238
- double sinSlope = sinZenithRad * std::sin ( slopeRad );
239
- double color0 = cosSlope + sinSlope * std::cos ( angle0Rad - aspectRad );
240
- double color1 = cosSlope + sinSlope * std::cos ( angle1Rad - aspectRad );
241
- double color2 = cosSlope + sinSlope * std::cos ( angle2Rad - aspectRad );
242
- double color3 = cosSlope + sinSlope * std::cos ( angle3Rad - aspectRad );
243
- grayValue = qBound ( 0.0 , 255 * ( weight0 * color0 + weight1 * color1 + weight2 * color2 + weight3 * color3 ) * 0.5 , 255.0 );
246
+ // Fast formula from GDAL DEM
247
+ const float xx = derX * derX;
248
+ const float yy = derY * derY;
249
+ const float xx_plus_yy = xx + yy;
250
+ // Flat?
251
+ if ( xx_plus_yy == 0.0 )
252
+ {
253
+ grayValue = qBound ( 0 .0f , static_cast <float >( 1.0 + sinZenithRad_mul_254 ), 255 .0f );
254
+ }
255
+ else
256
+ {
257
+ // ... then the shade value from different azimuth
258
+ float val225_mul_127 = sinZenithRad_mul_127 +
259
+ ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
260
+ val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127;
261
+ float val270_mul_127 = sinZenithRad_mul_127 -
262
+ derX * cos_alt_mul_z_mul_127;
263
+ val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127;
264
+ float val315_mul_127 = sinZenithRad_mul_127 +
265
+ ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127;
266
+ val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127;
267
+ float val360_mul_127 = sinZenithRad_mul_127 -
268
+ derY * cos_alt_mul_z_mul_127;
269
+ val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127;
270
+
271
+ // ... then the weighted shading
272
+ const float weight_225 = 0.5 * xx_plus_yy - derX * derY;
273
+ const float weight_270 = xx;
274
+ const float weight_315 = xx_plus_yy - weight_225;
275
+ const float weight_360 = yy;
276
+ const float cang_mul_127 = (
277
+ ( weight_225 * val225_mul_127 +
278
+ weight_270 * val270_mul_127 +
279
+ weight_315 * val315_mul_127 +
280
+ weight_360 * val360_mul_127 ) / xx_plus_yy ) /
281
+ ( 1 + square_z * xx_plus_yy );
282
+
283
+ grayValue = qBound ( 0 .0f , 1 .0f + cang_mul_127, 255 .0f );
284
+ }
244
285
}
245
286
246
- double currentAlpha = mOpacity ;
287
+ float currentAlpha = mOpacity ;
247
288
if ( mRasterTransparency )
248
289
{
249
290
currentAlpha = mRasterTransparency ->alphaValue ( x22, mOpacity * 255 ) / 255.0 ;
@@ -264,6 +305,17 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
264
305
}
265
306
}
266
307
308
+ #ifdef QGISDEBUG
309
+ if ( QgsSettings ().value ( QStringLiteral ( " Map/logCanvasRefreshEvent" ), false ).toBool () )
310
+ {
311
+ QgsMessageLog::logMessage ( QStringLiteral ( " CPU processing time for hillshade (%2 x %3 ): %4 ms" )
312
+ .arg ( width )
313
+ .arg ( height )
314
+ .arg ( std::chrono::duration_cast<std::chrono::milliseconds>( std::chrono::system_clock::now () - startTime ).count () ),
315
+ tr ( " Rendering" ) );
316
+ }
317
+ #endif
318
+
267
319
return outputBlock.release ();
268
320
}
269
321
0 commit comments