QGIS API Documentation  2.0.1-Dufour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgshillshadefilter.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgshillshadefilter.h - description
3  --------------------------------
4  begin : September 26th, 2011
5  copyright : (C) 2011 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
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 "qgshillshadefilter.h"
19 
20 QgsHillshadeFilter::QgsHillshadeFilter( const QString& inputFile, const QString& outputFile, const QString& outputFormat, double lightAzimuth,
21  double lightAngle )
22  : QgsDerivativeFilter( inputFile, outputFile, outputFormat )
23  , mLightAzimuth( lightAzimuth )
24  , mLightAngle( lightAngle )
25 {
26 }
27 
29 {
30 }
31 
32 float QgsHillshadeFilter::processNineCellWindow( float* x11, float* x21, float* x31,
33  float* x12, float* x22, float* x32,
34  float* x13, float* x23, float* x33 )
35 {
36  float derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33 );
37  float derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33 );
38 
39  if ( derX == mOutputNodataValue || derY == mOutputNodataValue )
40  {
41  return mOutputNodataValue;
42  }
43 
44  float zenith_rad = mLightAngle * M_PI / 180.0;
45  float slope_rad = atan( sqrt( derX * derX + derY * derY ) );
46  float azimuth_rad = mLightAzimuth * M_PI / 180.0;
47  float aspect_rad = 0;
48  if ( derX == 0 && derY == 0 ) //aspect undefined, take a neutral value. Better solutions?
49  {
50  aspect_rad = azimuth_rad / 2.0;
51  }
52  else
53  {
54  aspect_rad = M_PI + atan2( derX, derY );
55  }
56  return qMax( 0.0, 255.0 * (( cos( zenith_rad ) * cos( slope_rad ) ) + ( sin( zenith_rad ) * sin( slope_rad ) * cos( azimuth_rad - aspect_rad ) ) ) );
57 }