[opencl] Fix small OpenCL alg issues

From comparison tests with CPU results

+ some minor speed improvements

Author: elpaso

python/analysis/auto_generated/raster/qgsaspectfilter.sip.in

@@ -30,8 +30,6 @@ nodata value if not present or outside of the border. Must be implemented by sub
 %End
 
 
-
-
 };
 
 /************************************************************************

python/analysis/auto_generated/raster/qgsslopefilter.sip.in

@@ -29,6 +29,7 @@ Calculates output value from nine input values. The input values and the output
 nodata value if not present or outside of the border. Must be implemented by subclasses*
 %End
 
+
 };
 
 /************************************************************************

resources/opencl_programs/hillshade.cl

@@ -1,8 +1,11 @@
 #include "calcfirstder.cl"
+<<<<<<< 30a62e1addd2541bdbecda4398381772eb8034d4
 
 #ifndef M_PI
     #define M_PI 3.14159265358979323846
 #endif
+=======
+>>>>>>> [opencl] Fix small OpenCL alg issues
 
 __kernel void processNineCellWindow( __global float *scanLine1,
                                      __global float *scanLine2,
@@ -39,6 +42,7 @@ __kernel void processNineCellWindow( __global float *scanLine1,
       else
       {
 
+<<<<<<< 30a62e1addd2541bdbecda4398381772eb8034d4
         float slope_rad = sqrt( derX * derX + derY * derY );
         slope_rad = atan( slope_rad );
         float aspect_rad;
@@ -50,6 +54,23 @@ __kernel void processNineCellWindow( __global float *scanLine1,
         resultLine[i] = max(0.0f, 255.0f * ( ( rasterParams[5] * cos( slope_rad ) ) +
                                             (  rasterParams[6] * sin( slope_rad ) *
                                               cos( rasterParams[7] - aspect_rad ) ) ) );
+=======
+        float zenith_rad = rasterParams[5];
+        float azimuth_rad = rasterParams[6];
+        float slope_rad = sqrt( derX * derX + derY * derY );
+
+
+        slope_rad = atan( slope_rad );
+        float aspect_rad;
+        if ( derX == 0.0f && derY == 0.0f)
+          aspect_rad = azimuth_rad / 2.0f;
+        else
+          aspect_rad = M_PI + atan2( derX, derY );
+
+        resultLine[i] = max(0.0f, 255.0f * ( ( cos( zenith_rad ) * cos( slope_rad ) ) +
+                                            ( sin( zenith_rad ) * sin( slope_rad ) *
+                                              cos( azimuth_rad - aspect_rad ) ) ) );
+>>>>>>> [opencl] Fix small OpenCL alg issues
 
       }
     }

src/analysis/raster/qgsalignraster.cpp

@@ -615,3 +615,4 @@ double QgsAlignRaster::RasterInfo::identify( double mx, double my )
 
   return value;
 }
+

src/analysis/raster/qgsaspectfilter.h

@@ -37,14 +37,15 @@ class ANALYSIS_EXPORT QgsAspectFilter: public QgsDerivativeFilter
                                  float *x13, float *x23, float *x33 ) override;
 
 
-
-
-    // QgsNineCellFilter interface
+#ifdef HAVE_OPENCL
   private:
+
     const QString openClProgramBaseName() const override
     {
       return QStringLiteral( "aspect" );
     }
+#endif
+
 };
 
 #endif // QGSASPECTFILTER_H

src/analysis/raster/qgshillshadefilter.cpp

@@ -21,15 +21,19 @@
 QgsHillshadeFilter::QgsHillshadeFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat, double lightAzimuth,
                                         double lightAngle )
   : QgsDerivativeFilter( inputFile, outputFile, outputFormat )
-  , mLightAzimuth( lightAzimuth )
-  , mLightAngle( lightAngle )
+  , mLightAzimuth( static_cast<float>( lightAzimuth ) )
+  , mLightAngle( static_cast<float>( lightAngle ) )
+  , mCosZenithRad( std::cos( static_cast<float>( lightAngle * M_PI ) / 180.0f ) )
+  , mSinZenithRad( std::sin( static_cast<float>( lightAngle * M_PI ) / 180.0f ) )
+  , mAzimuthRad( static_cast<float>( lightAzimuth * M_PI ) / 180.0f )
 {
 }
 
 float QgsHillshadeFilter::processNineCellWindow( float *x11, float *x21, float *x31,
     float *x12, float *x22, float *x32,
     float *x13, float *x23, float *x33 )
 {
+
   float derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33 );
   float derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33 );
 
@@ -38,17 +42,43 @@ float QgsHillshadeFilter::processNineCellWindow( float *x11, float *x21, float *
     return mOutputNodataValue;
   }
 
-  float zenith_rad = mLightAngle * M_PI / 180.0;
   float slope_rad = std::atan( std::sqrt( derX * derX + derY * derY ) );
-  float azimuth_rad = mLightAzimuth * M_PI / 180.0;
   float aspect_rad = 0;
   if ( derX == 0 && derY == 0 ) //aspect undefined, take a neutral value. Better solutions?
   {
-    aspect_rad = azimuth_rad / 2.0;
+    aspect_rad = mAzimuthRad / 2.0f;
   }
   else
   {
     aspect_rad = M_PI + std::atan2( derX, derY );
   }
-  return std::max( 0.0, 255.0 * ( ( std::cos( zenith_rad ) * std::cos( slope_rad ) ) + ( std::sin( zenith_rad ) * std::sin( slope_rad ) * std::cos( azimuth_rad - aspect_rad ) ) ) );
+  return std::max( 0.0f, 255.0f * ( ( mCosZenithRad * std::cos( slope_rad ) ) +
+                                    ( mSinZenithRad * std::sin( slope_rad ) *
+                                      std::cos( mAzimuthRad - aspect_rad ) ) ) );
+}
+
+void QgsHillshadeFilter::setLightAzimuth( float azimuth )
+{
+  mLightAzimuth = azimuth;
+  mAzimuthRad = azimuth * static_cast<float>( M_PI ) / 180.0f;
+}
+
+void QgsHillshadeFilter::setLightAngle( float angle )
+{
+  mLightAngle = angle;
+  mCosZenithRad = std::cos( angle * static_cast<float>( M_PI ) / 180.0f );
+  mSinZenithRad = std::sin( angle * static_cast<float>( M_PI ) / 180.0f );
 }
+
+#ifdef HAVE_OPENCL
+
+void QgsHillshadeFilter::addExtraRasterParams( std::vector<float> &params )
+{
+
+  params.push_back( mCosZenithRad ); // cos_zenith_rad 5
+  params.push_back( mSinZenithRad ); // sin_zenith_rad 6
+  params.push_back( mAzimuthRad ); // azimuth_rad 7
+
+}
+
+#endif

src/analysis/raster/qgshillshadefilter.h

@@ -39,13 +39,33 @@ class ANALYSIS_EXPORT QgsHillshadeFilter: public QgsDerivativeFilter
                                  float *x13, float *x23, float *x33 ) override;
 
     float lightAzimuth() const { return mLightAzimuth; }
-    void setLightAzimuth( float azimuth ) { mLightAzimuth = azimuth; }
+    void setLightAzimuth( float azimuth );
     float lightAngle() const { return mLightAngle; }
-    void setLightAngle( float angle ) { mLightAngle = angle; }
+    void setLightAngle( float angle );
 
   private:
+
+#ifdef HAVE_OPENCL
+
+    const QString openClProgramBaseName() const override
+    {
+      return QStringLiteral( "hillshade" );
+    }
+#endif
+
     float mLightAzimuth;
     float mLightAngle;
+    // Precalculate for speed:
+    float mCosZenithRad;
+    float mSinZenithRad;
+    float mAzimuthRad;
+
+
+#ifdef HAVE_OPENCL
+  private:
+
+    void addExtraRasterParams( std::vector<float> &params ) override;
+#endif
 
 };
 

src/analysis/raster/qgsninecellfilter.cpp

@@ -21,12 +21,18 @@
 #include "qgsfeedback.h"
 #include "qgsogrutils.h"
 #include "qgsmessagelog.h"
+
+#ifdef HAVE_OPENCL
+#include "qgsopenclutils.h"
+#endif
+
 #include <QFile>
 #include <QDebug>
 #include <QFileInfo>
 #include <iterator>
 
 
+
 QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
   : mInputFile( inputFile )
   , mOutputFile( outputFile )
@@ -157,6 +163,7 @@ gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDa
   return outputDataset;
 }
 
+#ifdef HAVE_OPENCL
 
 int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *feedback )
 {
@@ -296,7 +303,8 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
         }
         queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
       }
-      else // Overwrite from input, skip first and last
+      else // Read line i + 1 and put it into scanline 3
+        // Overwrite from input, skip first and last
       {
         if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
         {
@@ -334,7 +342,7 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
   }
   return 0;
 }
-
+#endif
 
 int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
 {
@@ -393,7 +401,7 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
   float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
 
   //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
-  for ( int i = 0; i < ySize; ++i )
+  for ( int yIndex = 0; yIndex < ySize; ++yIndex )
   {
     if ( feedback && feedback->isCanceled() )
     {
@@ -402,10 +410,10 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
 
     if ( feedback )
     {
-      feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize );
+      feedback->setProgress( 100.0 * static_cast< double >( yIndex ) / ySize );
     }
 
-    if ( i == 0 )
+    if ( yIndex == 0 )
     {
       //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
       for ( int a = 0; a < xSize + 2 ; ++a )
@@ -428,7 +436,7 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
     }
 
     // Read scanline 3
-    if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
+    if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values
     {
       for ( int a = 0; a < xSize + 2; ++a )
       {
@@ -437,7 +445,7 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
     }
     else
     {
-      if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
+      if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
       {
         QgsDebugMsg( "Raster IO Error" );
       }
@@ -450,15 +458,16 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
 
 
     // j is the x axis index, skip 0 and last cell that have been filled with nodata
-    for ( int j = 0; j < xSize ; ++j )
+    for ( int xIndex = 0; xIndex < xSize ; ++xIndex )
     {
-      resultLine[ j ] = processNineCellWindow( &scanLine1[ j ], &scanLine1[ j + 1 ], &scanLine1[ j + 2 ],
-                        &scanLine2[ j ], &scanLine2[ j + 1 ], &scanLine2[ j + 2 ],
-                        &scanLine3[ j ], &scanLine3[ j + 1 ], &scanLine3[ j + 2 ] );
+      // cells(x, y) x11, x21, x31, x12, x22, x32, x13, x23, x33
+      resultLine[ xIndex ] = processNineCellWindow( &scanLine1[ xIndex ], &scanLine1[ xIndex + 1 ], &scanLine1[ xIndex + 2 ],
+                             &scanLine2[ xIndex ], &scanLine2[ xIndex + 1 ], &scanLine2[ xIndex + 2 ],
+                             &scanLine3[ xIndex ], &scanLine3[ xIndex + 1 ], &scanLine3[ xIndex + 2 ] );
 
     }
 
-    if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
+    if ( GDALRasterIO( outputRasterBand, GF_Write, 0, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
     {
       QgsDebugMsg( "Raster IO Error" );
     }

src/analysis/raster/qgsninecellfilter.h

@@ -22,7 +22,6 @@
 #include "gdal.h"
 #include "qgis_analysis.h"
 #include "qgsogrutils.h"
-#include "qgsopenclutils.h"
 
 class QgsFeedback;
 

src/analysis/raster/qgsruggednessfilter.cpp

@@ -24,11 +24,6 @@ QgsRuggednessFilter::QgsRuggednessFilter( const QString &inputFile, const QStrin
 
 }
 
-QgsRuggednessFilter::QgsRuggednessFilter()
-  : QgsNineCellFilter( QString(), QString(), QString() )
-{
-
-}
 
 float QgsRuggednessFilter::processNineCellWindow( float *x11, float *x21, float *x31,
     float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 )