Skip to content

Commit 1c9b65d

Browse files
authoredJul 19, 2018
Merge pull request #7411 from elpaso/hillshade-fast
Hillshade renderer speed improvements
2 parents e6af079 + 3138339 commit 1c9b65d

File tree

1 file changed

+100
-48
lines changed

1 file changed

+100
-48
lines changed
 

‎src/core/raster/qgshillshaderenderer.cpp

Lines changed: 100 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -22,8 +22,14 @@
2222
#include "qgsrasterinterface.h"
2323
#include "qgsrasterblock.h"
2424
#include "qgsrectangle.h"
25+
#include "qgssettings.h"
2526
#include <memory>
2627

28+
#ifdef QGISDEBUG
29+
#include "qgsmessagelog.h"
30+
#include <chrono>
31+
#endif
32+
2733
QgsHillshadeRenderer::QgsHillshadeRenderer( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ):
2834
QgsRasterRenderer( input, QStringLiteral( "hillshade" ) )
2935
, mBand( band )
@@ -122,21 +128,34 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
122128
return outputBlock.release();
123129
}
124130

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;
137152

138153
QRgb myDefaultColor = NODATA_COLOR;
139154

155+
#ifdef QGISDEBUG
156+
std::chrono::time_point<std::chrono::system_clock> startTime( std::chrono::system_clock::now() );
157+
#endif
158+
140159
for ( qgssize i = 0; i < ( qgssize )height; i++ )
141160
{
142161

@@ -182,15 +201,15 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
182201
jRight = j;
183202
}
184203

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;
194213

195214
// This is center cell. It is not nodata. Use this in place of nodata neighbors
196215
x22 = inputBlock->value( i, j );
@@ -207,43 +226,65 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
207226
x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight );
208227
x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight );
209228

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 );
215231

216-
217-
double grayValue;
232+
float grayValue;
218233
if ( !mMultiDirectional )
219234
{
220235
// 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 );
224242
}
225243
else
226244
{
227245
// 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+
}
244285
}
245286

246-
double currentAlpha = mOpacity;
287+
float currentAlpha = mOpacity;
247288
if ( mRasterTransparency )
248289
{
249290
currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0;
@@ -264,6 +305,17 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext
264305
}
265306
}
266307

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+
267319
return outputBlock.release();
268320
}
269321

0 commit comments

Comments
 (0)
Please sign in to comment.