QGIS API Documentation  3.0.2-Girona (307d082)
qgshillshaderenderer.cpp
Go to the documentation of this file.
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 <memory>
26 
27 QgsHillshadeRenderer::QgsHillshadeRenderer( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ):
28  QgsRasterRenderer( input, QStringLiteral( "hillshade" ) )
29  , mBand( band )
30  , mZFactor( 1 )
31  , mLightAngle( lightAngle )
32  , mLightAzimuth( lightAzimuth )
33  , mMultiDirectional( false )
34 {
35 
36 }
37 
39 {
40  QgsHillshadeRenderer *r = new QgsHillshadeRenderer( nullptr, mBand, mLightAzimuth, mLightAngle );
41  r->copyCommonProperties( this );
42 
43  r->setZFactor( mZFactor );
44  r->setMultiDirectional( mMultiDirectional );
45  return r;
46 }
47 
49 {
50  if ( elem.isNull() )
51  {
52  return nullptr;
53  }
54 
55  int band = elem.attribute( QStringLiteral( "band" ), QStringLiteral( "0" ) ).toInt();
56  double azimuth = elem.attribute( QStringLiteral( "azimuth" ), QStringLiteral( "315" ) ).toDouble();
57  double angle = elem.attribute( QStringLiteral( "angle" ), QStringLiteral( "45" ) ).toDouble();
58  double zFactor = elem.attribute( QStringLiteral( "zfactor" ), QStringLiteral( "1" ) ).toDouble();
59  bool multiDirectional = elem.attribute( QStringLiteral( "multidirection" ), QStringLiteral( "0" ) ).toInt();
60  QgsHillshadeRenderer *r = new QgsHillshadeRenderer( input, band, azimuth, angle );
61  r->readXml( elem );
62 
63  r->setZFactor( zFactor );
64  r->setMultiDirectional( multiDirectional );
65  return r;
66 }
67 
68 void QgsHillshadeRenderer::writeXml( QDomDocument &doc, QDomElement &parentElem ) const
69 {
70  if ( parentElem.isNull() )
71  {
72  return;
73  }
74 
75  QDomElement rasterRendererElem = doc.createElement( QStringLiteral( "rasterrenderer" ) );
76  _writeXml( doc, rasterRendererElem );
77 
78  rasterRendererElem.setAttribute( QStringLiteral( "band" ), mBand );
79  rasterRendererElem.setAttribute( QStringLiteral( "azimuth" ), QString::number( mLightAzimuth ) );
80  rasterRendererElem.setAttribute( QStringLiteral( "angle" ), QString::number( mLightAngle ) );
81  rasterRendererElem.setAttribute( QStringLiteral( "zfactor" ), QString::number( mZFactor ) );
82  rasterRendererElem.setAttribute( QStringLiteral( "multidirection" ), QString::number( mMultiDirectional ) );
83  parentElem.appendChild( rasterRendererElem );
84 }
85 
86 QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
87 {
88  Q_UNUSED( bandNo );
89  std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock() );
90  if ( !mInput )
91  {
92  QgsDebugMsg( "No input raster!" );
93  return outputBlock.release();
94  }
95 
96  std::shared_ptr< QgsRasterBlock > inputBlock( mInput->block( mBand, extent, width, height, feedback ) );
97 
98  if ( !inputBlock || inputBlock->isEmpty() )
99  {
100  QgsDebugMsg( "No raster data!" );
101  return outputBlock.release();
102  }
103 
104  std::shared_ptr< QgsRasterBlock > alphaBlock;
105 
106  if ( mAlphaBand > 0 && mBand != mAlphaBand )
107  {
108  alphaBlock.reset( mInput->block( mAlphaBand, extent, width, height, feedback ) );
109  if ( !alphaBlock || alphaBlock->isEmpty() )
110  {
111  // TODO: better to render without alpha
112  return outputBlock.release();
113  }
114  }
115  else if ( mAlphaBand > 0 )
116  {
117  alphaBlock = inputBlock;
118  }
119 
120  if ( !outputBlock->reset( Qgis::ARGB32_Premultiplied, width, height ) )
121  {
122  return outputBlock.release();
123  }
124 
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;
137 
138  QRgb myDefaultColor = NODATA_COLOR;
139 
140  for ( qgssize i = 0; i < ( qgssize )height; i++ )
141  {
142 
143  for ( qgssize j = 0; j < ( qgssize )width; j++ )
144  {
145 
146  if ( inputBlock->isNoData( i, j ) )
147  {
148  outputBlock->setColor( i, j, myDefaultColor );
149  continue;
150  }
151 
152  qgssize iUp, iDown, jLeft, jRight;
153  if ( i == 0 )
154  {
155  iUp = i;
156  iDown = i + 1;
157  }
158  else if ( i < ( qgssize )height - 1 )
159  {
160  iUp = i - 1;
161  iDown = i + 1;
162  }
163  else
164  {
165  iUp = i - 1;
166  iDown = i;
167  }
168 
169  if ( j == 0 )
170  {
171  jLeft = j;
172  jRight = j + 1;
173  }
174  else if ( j < ( qgssize )width - 1 )
175  {
176  jLeft = j - 1;
177  jRight = j + 1;
178  }
179  else
180  {
181  jLeft = j - 1;
182  jRight = j;
183  }
184 
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;
194 
195  // This is center cell. It is not nodata. Use this in place of nodata neighbors
196  x22 = inputBlock->value( i, j );
197 
198  x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft );
199  x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft );
200  x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft );
201 
202  x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j );
203  // x22
204  x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j );
205 
206  x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight );
207  x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight );
208  x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight );
209 
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 );
215 
216 
217  double grayValue;
218  if ( !mMultiDirectional )
219  {
220  // 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 );
224  }
225  else
226  {
227  // 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 );
244  }
245 
246  double currentAlpha = mOpacity;
247  if ( mRasterTransparency )
248  {
249  currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0;
250  }
251  if ( mAlphaBand > 0 )
252  {
253  currentAlpha *= alphaBlock->value( i ) / 255.0;
254  }
255 
256  if ( qgsDoubleNear( currentAlpha, 1.0 ) )
257  {
258  outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) );
259  }
260  else
261  {
262  outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
263  }
264  }
265  }
266 
267  return outputBlock.release();
268 }
269 
271 {
272  QList<int> bandList;
273  if ( mBand != -1 )
274  {
275  bandList << mBand;
276  }
277  return bandList;
278 
279 }
280 
282 {
283  if ( bandNo > mInput->bandCount() || bandNo <= 0 )
284  {
285  return;
286  }
287  mBand = bandNo;
288 }
289 
290 double QgsHillshadeRenderer::calcFirstDerX( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
291 {
292  Q_UNUSED( x12 );
293  Q_UNUSED( x22 );
294  Q_UNUSED( x32 );
295  return ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellsize );
296 }
297 
298 double QgsHillshadeRenderer::calcFirstDerY( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
299 {
300  Q_UNUSED( x21 );
301  Q_UNUSED( x22 );
302  Q_UNUSED( x23 );
303  return ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellsize );
304 }
305 
306 
307 
308 
virtual int bandCount() const =0
Get number of bands.
int alphaValue(double, int globalTransparency=255) const
Returns the transparency value for a single value Pixel.
A rectangle specified with double values.
Definition: qgsrectangle.h:39
void setMultiDirectional(bool isMultiDirectional)
Sets whether to render using a multi-directional hillshade algorithm.
virtual QgsRectangle extent() const
Get the extent of the interface.
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
virtual QgsRasterInterface * input() const
Current input.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:251
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:786
void setBand(int bandNo)
Sets the band used by the renderer.
QgsRasterTransparency * mRasterTransparency
Raster transparency per color or value. Overwrites global alpha value.
void setZFactor(double zfactor)
Set the Z scaling factor of the result image.
int band() const
Returns the band used by the renderer.
static const QRgb NODATA_COLOR
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:106
QList< int > usesBands() const override
Returns a list of band numbers used by the renderer.
Raster data container.
bool multiDirectional() const
Returns true if the renderer is using multi-directional hillshading.
virtual QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr)=0
Read block of data using given extent and size.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:142
void _writeXml(QDomDocument &doc, QDomElement &rasterRendererElem) const
Write upper class info into rasterrenderer element (called by writeXml method of subclasses) ...
void copyCommonProperties(const QgsRasterRenderer *other, bool copyMinMaxOrigin=true)
Copies common properties like opacity / transparency data from other renderer.
A renderer for generating live hillshade models.
QgsHillshadeRenderer(QgsRasterInterface *input, int band, double lightAzimuth, double lightAltitude)
A renderer for generating live hillshade models.
int mAlphaBand
Read alpha value from band.
void readXml(const QDomElement &rendererElem) override
Sets base class members from xml. Usually called from create() methods of subclasses.
Base class for processing filters like renderers, reprojector, resampler etc.
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:488
QgsHillshadeRenderer * clone() const override
Clone itself, create deep copy.
void writeXml(QDomDocument &doc, QDomElement &parentElem) const override
Write base class members to xml.
double mOpacity
Global alpha value (0-1)
double zFactor() const
Returns the Z scaling factor.
double azimuth() const
Returns the direction of the light over the raster between 0-360.
static QgsRasterRenderer * create(const QDomElement &elem, QgsRasterInterface *input)
Factory method to create a new renderer.
QgsRasterInterface * mInput
Feedback object tailored for raster block reading.
Raster renderer pipe that applies colors to a raster.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:149