Quantum GIS API Documentation
1.8
|
00001 /*************************************************************************** 00002 qgshillshadefilter.h - description 00003 -------------------------------- 00004 begin : September 26th, 2011 00005 copyright : (C) 2011 by Marco Hugentobler 00006 email : marco dot hugentobler at sourcepole dot ch 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 ***************************************************************************/ 00017 00018 #include "qgshillshadefilter.h" 00019 00020 QgsHillshadeFilter::QgsHillshadeFilter( const QString& inputFile, const QString& outputFile, const QString& outputFormat, double lightAzimuth, 00021 double lightAngle ) 00022 : QgsDerivativeFilter( inputFile, outputFile, outputFormat ) 00023 , mLightAzimuth( lightAzimuth ) 00024 , mLightAngle( lightAngle ) 00025 { 00026 } 00027 00028 QgsHillshadeFilter::~QgsHillshadeFilter() 00029 { 00030 } 00031 00032 float QgsHillshadeFilter::processNineCellWindow( float* x11, float* x21, float* x31, 00033 float* x12, float* x22, float* x32, 00034 float* x13, float* x23, float* x33 ) 00035 { 00036 float derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33 ); 00037 float derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33 ); 00038 00039 if ( derX == mOutputNodataValue || derY == mOutputNodataValue ) 00040 { 00041 return mOutputNodataValue; 00042 } 00043 00044 float zenith_rad = mLightAngle * M_PI / 180.0; 00045 float slope_rad = atan( sqrt( derX * derX + derY * derY ) ); 00046 float azimuth_rad = mLightAzimuth * M_PI / 180.0; 00047 float aspect_rad = 0; 00048 if ( derX == 0 && derY == 0 ) //aspect undefined, take a neutral value. Better solutions? 00049 { 00050 aspect_rad = azimuth_rad / 2.0; 00051 } 00052 else 00053 { 00054 aspect_rad = M_PI + atan2( derX, derY ); 00055 } 00056 return qMax( 0.0, 255.0 * (( cos( zenith_rad ) * cos( slope_rad ) ) + ( sin( zenith_rad ) * sin( slope_rad ) * cos( azimuth_rad - aspect_rad ) ) ) ); 00057 }