QGIS API Documentation  2.18.21-Las Palmas (9fba24a)
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 
26 
27 QgsHillshadeRenderer::QgsHillshadeRenderer( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ):
28  QgsRasterRenderer( input, "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( "band", "0" ).toInt();
56  double azimuth = elem.attribute( "azimuth", "315" ).toDouble();
57  double angle = elem.attribute( "angle", "45" ).toDouble();
58  double zFactor = elem.attribute( "zfactor", "1" ).toDouble();
59  bool multiDirectional = elem.attribute( "multidirection", "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 
69 {
70  if ( parentElem.isNull() )
71  {
72  return;
73  }
74 
75  QDomElement rasterRendererElem = doc.createElement( "rasterrenderer" );
76  _writeXML( doc, rasterRendererElem );
77 
78  rasterRendererElem.setAttribute( "band", mBand );
79  rasterRendererElem.setAttribute( "azimuth", QString::number( mLightAzimuth ) );
80  rasterRendererElem.setAttribute( "angle", QString::number( mLightAngle ) );
81  rasterRendererElem.setAttribute( "zfactor", QString::number( mZFactor ) );
82  rasterRendererElem.setAttribute( "multidirection", QString::number( mMultiDirectional ) );
83  parentElem.appendChild( rasterRendererElem );
84 }
85 
86 QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &extent, int width, int height )
87 {
88  return block2( bandNo, extent, width, height );
89 }
90 
91 QgsRasterBlock *QgsHillshadeRenderer::block2( int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback* feedback )
92 {
93  Q_UNUSED( bandNo );
94  QgsRasterBlock *outputBlock = new QgsRasterBlock();
95  if ( !mInput )
96  {
97  QgsDebugMsg( "No input raster!" );
98  return outputBlock;
99  }
100 
101  QgsRasterBlock *inputBlock = mInput->block2( mBand, extent, width, height, feedback );
102 
103  if ( !inputBlock || inputBlock->isEmpty() )
104  {
105  QgsDebugMsg( "No raster data!" );
106  delete inputBlock;
107  return outputBlock;
108  }
109 
110  QgsRasterBlock *alphaBlock = nullptr;
111 
112  if ( mAlphaBand > 0 && mBand != mAlphaBand )
113  {
114 
115  alphaBlock = mInput->block2( mAlphaBand, extent, width, height, feedback );
116  if ( !alphaBlock || alphaBlock->isEmpty() )
117  {
118  // TODO: better to render without alpha
119  delete inputBlock;
120  delete alphaBlock;
121  return outputBlock;
122  }
123  }
124  else if ( mAlphaBand > 0 )
125  {
126  alphaBlock = inputBlock;
127  }
128 
129  if ( !outputBlock->reset( QGis::ARGB32_Premultiplied, width, height ) )
130  {
131  delete inputBlock;
132  delete alphaBlock;
133  return outputBlock;
134  }
135 
136 
137 
138  double cellXSize = extent.width() / double( width );
139  double cellYSize = extent.height() / double( height );
140  double zenithRad = qMax( 0.0, 90 - mLightAngle ) * M_PI / 180.0;
141  double azimuthRad = -1 * mLightAzimuth * M_PI / 180.0;
142  double cosZenithRad = cos( zenithRad );
143  double sinZenithRad = sin( zenithRad );
144 
145  // Multi direction hillshade: http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
146  double angle0Rad = ( -1 * mLightAzimuth - 45 - 45 * 0.5 ) * M_PI / 180.0;
147  double angle1Rad = ( -1 * mLightAzimuth - 45 * 0.5 ) * M_PI / 180.0;
148  double angle2Rad = ( -1 * mLightAzimuth + 45 * 0.5 ) * M_PI / 180.0;
149  double angle3Rad = ( -1 * mLightAzimuth + 45 + 45 * 0.5 ) * M_PI / 180.0;
150 
151  QRgb myDefaultColor = NODATA_COLOR;
152 
153  for ( qgssize i = 0; i < ( qgssize )height; i++ )
154  {
155 
156  for ( qgssize j = 0; j < ( qgssize )width; j++ )
157  {
158 
159  if ( inputBlock->isNoData( i, j ) )
160  {
161  outputBlock->setColor( i, j, myDefaultColor );
162  continue;
163  }
164 
165  qgssize iUp, iDown, jLeft, jRight;
166  if ( i == 0 )
167  {
168  iUp = i;
169  iDown = i + 1;
170  }
171  else if ( i < ( qgssize )height - 1 )
172  {
173  iUp = i - 1;
174  iDown = i + 1;
175  }
176  else
177  {
178  iUp = i - 1;
179  iDown = i;
180  }
181 
182  if ( j == 0 )
183  {
184  jLeft = j;
185  jRight = j + 1;
186  }
187  else if ( j < ( qgssize )width - 1 )
188  {
189  jLeft = j - 1;
190  jRight = j + 1;
191  }
192  else
193  {
194  jLeft = j - 1;
195  jRight = j;
196  }
197 
198  double x11;
199  double x21;
200  double x31;
201  double x12;
202  double x22; // Working cell
203  double x32;
204  double x13;
205  double x23;
206  double x33;
207 
208  // This is center cell. It is not nodata. Use this in place of nodata neighbors
209  x22 = inputBlock->value( i, j );
210 
211  x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft );
212  x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft );
213  x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft );
214 
215  x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j );
216  // x22
217  x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j );
218 
219  x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight );
220  x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight );
221  x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight );
222 
223  double derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize );
224  double derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize );
225 
226  double slopeRad = atan( mZFactor * sqrt( derX * derX + derY * derY ) );
227  double aspectRad = atan2( derX, -derY );
228 
229 
230  double grayValue;
231  if ( !mMultiDirectional )
232  {
233  // Standard single direction hillshade
234  grayValue = qBound( 0.0, 255.0 * ( cosZenithRad * cos( slopeRad )
235  + sinZenithRad * sin( slopeRad )
236  * cos( azimuthRad - aspectRad ) ) , 255.0 );
237  }
238  else
239  {
240  // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
241  double weight0 = sin( aspectRad - angle0Rad );
242  double weight1 = sin( aspectRad - angle1Rad );
243  double weight2 = sin( aspectRad - angle2Rad );
244  double weight3 = sin( aspectRad - angle3Rad );
245  weight0 = weight0 * weight0;
246  weight1 = weight1 * weight1;
247  weight2 = weight2 * weight2;
248  weight3 = weight3 * weight3;
249 
250  double cosSlope = cosZenithRad * cos( slopeRad );
251  double sinSlope = sinZenithRad * sin( slopeRad );
252  double color0 = cosSlope + sinSlope * cos( angle0Rad - aspectRad ) ;
253  double color1 = cosSlope + sinSlope * cos( angle1Rad - aspectRad ) ;
254  double color2 = cosSlope + sinSlope * cos( angle2Rad - aspectRad ) ;
255  double color3 = cosSlope + sinSlope * cos( angle3Rad - aspectRad ) ;
256  grayValue = qBound( 0.0, 255 * ( weight0 * color0 + weight1 * color1 + weight2 * color2 + weight3 * color3 ) * 0.5, 255.0 );
257  }
258 
259  double currentAlpha = mOpacity;
260  if ( mRasterTransparency )
261  {
262  currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0;
263  }
264  if ( mAlphaBand > 0 )
265  {
266  currentAlpha *= alphaBlock->value( i ) / 255.0;
267  }
268 
269  if ( qgsDoubleNear( currentAlpha, 1.0 ) )
270  {
271  outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) );
272  }
273  else
274  {
275  outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) );
276  }
277  }
278  }
279 
280  delete inputBlock;
281  if ( mAlphaBand > 0 && mBand != mAlphaBand )
282  {
283  delete alphaBlock;
284  }
285  return outputBlock;
286 }
287 
289 {
290  QList<int> bandList;
291  if ( mBand != -1 )
292  {
293  bandList << mBand;
294  }
295  return bandList;
296 
297 }
298 
300 {
301  if ( bandNo > mInput->bandCount() || bandNo <= 0 )
302  {
303  return;
304  }
305  mBand = bandNo;
306 }
307 
308 double QgsHillshadeRenderer::calcFirstDerX( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
309 {
310  Q_UNUSED( x12 );
311  Q_UNUSED( x22 );
312  Q_UNUSED( x32 );
313  return (( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * cellsize );
314 }
315 
316 double QgsHillshadeRenderer::calcFirstDerY( double x11, double x21, double x31, double x12, double x22, double x32, double x13, double x23, double x33, double cellsize )
317 {
318  Q_UNUSED( x21 );
319  Q_UNUSED( x22 );
320  Q_UNUSED( x23 );
321  return (( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -cellsize );
322 }
323 
324 
325 
326 
virtual int bandCount() const =0
Get number of bands.
void writeXML(QDomDocument &doc, QDomElement &parentElem) const override
Write base class members to xml.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QgsRasterBlock * block2(int bandNo, QgsRectangle const &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
QDomNode appendChild(const QDomNode &newChild)
void setMultiDirectional(bool isMultiDirectional)
Sets whether to render using a multi-directional hillshade algorithm.
QString attribute(const QString &name, const QString &defValue) const
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
int alphaValue(double, int theGlobalTransparency=255) const
Returns the transparency value for a single value Pixel.
virtual QgsRasterInterface * input() const
Current input.
void copyCommonProperties(const QgsRasterRenderer *other)
Copies common properties like opacity / transparency data from other renderer.
double toDouble(bool *ok) const
bool isNoData(int row, int column)
Check if value at position is no data.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:353
void readXML(const QDomElement &rendererElem) override
Sets base class members from xml.
virtual QgsRasterBlock * block2(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr)
Read block of data using given extent and size.
void setBand(int bandNo)
Sets the band used by the renderer.
QgsRasterTransparency * mRasterTransparency
Raster transparency per color or value.
void setZFactor(double zfactor)
Set the Z scaling factor of the result image.
int band() const
Returns the band used by the renderer.
bool setColor(int row, int column, QRgb color)
Set color on position.
static const QRgb NODATA_COLOR
QString number(int n, int base)
QList< int > usesBands() const override
Returns a list of band numbers used by the renderer.
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:150
Raster data container.
bool multiDirectional() const
Returns true if the renderer is using multi-directional hillshading.
void setAttribute(const QString &name, const QString &value)
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
int toInt(bool *ok, int base) const
#define M_PI
A renderer for generating live hillshade models.
bool isEmpty() const
Returns true if block is empty, i.e.
void _writeXML(QDomDocument &doc, QDomElement &rasterRendererElem) const
Write upper class info into rasterrenderer element (called by writeXML method of subclasses) ...
QgsHillshadeRenderer(QgsRasterInterface *input, int band, double lightAzimuth, double lightAltitude)
A renderer for generating live hillshade models.
int mAlphaBand
Read alpha value from band.
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:500
double ANALYSIS_EXPORT angle(Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
QgsHillshadeRenderer * clone() const override
Clone itself, create deep copy.
bool reset(QGis::DataType theDataType, int theWidth, int theHeight)
Reset block.
bool isNull() const
virtual QgsRectangle extent()
Get the extent of the interface.
double value(int row, int column) const
Read a single value if type of block is numeric.
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.
QDomElement createElement(const QString &tagName)
static QgsRasterRenderer * create(const QDomElement &elem, QgsRasterInterface *input)
Factory method to create a new renderer.
QgsRasterInterface * mInput
Feedback object tailored for raster block reading.
QgsRasterBlock * block(int bandNo, QgsRectangle const &extent, int width, int height) override
Read block of data using given extent and size.
Raster renderer pipe that applies colors to a raster.
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212