Use OpenCL command queue

resources/CMakeLists.txt

@@ -11,6 +11,9 @@ INSTALL(DIRECTORY themes DESTINATION ${QGIS_DATA_DIR}/resources)
 INSTALL(DIRECTORY data DESTINATION ${QGIS_DATA_DIR}/resources)
 INSTALL(DIRECTORY metadata-ISO DESTINATION ${QGIS_DATA_DIR}/resources)
 INSTALL(DIRECTORY palettes DESTINATION ${QGIS_DATA_DIR}/resources)
+IF (HAVE_OPENCL)
+INSTALL(DIRECTORY opencl_programs DESTINATION ${QGIS_DATA_DIR}/resources)
+ENDIF (HAVE_OPENCL)
 
 IF (WITH_SERVER)
         INSTALL(DIRECTORY server DESTINATION ${QGIS_DATA_DIR}/resources)

resources/opencl_programs/aspect.cl

@@ -0,0 +1,45 @@
+#include "calcfirstder.cl"
+
+__kernel void processNineCellWindow( __global float *scanLine1,
+                                     __global float *scanLine2,
+                                     __global float *scanLine3,
+                                     __global float *resultLine,
+                                     __global float *rasterParams
+                       ) {
+
+  // Get the index of the current element
+  const int i = get_global_id(0);
+
+  // Do the operation
+  //return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * mCellSizeX))
+  float derX = calcFirstDer(   scanLine1[i],   scanLine2[i],   scanLine3[i],
+                               scanLine1[i+1], scanLine2[i+1], scanLine3[i+1],
+                               scanLine1[i+2], scanLine2[i+2], scanLine3[i+2],
+                               rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3]
+                             );
+  //return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * mCellSizeY));
+  float derY = calcFirstDer(   scanLine1[i+2], scanLine1[i+1], scanLine1[i],
+                               scanLine2[i+2], scanLine2[i+1], scanLine2[i],
+                               scanLine3[i+2], scanLine3[i+1], scanLine3[i],
+                               rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]
+                             );
+
+
+  if ( derX == rasterParams[1] || derY == rasterParams[1] ||
+       ( derX == 0.0f && derY == 0.0f) )
+  {
+    resultLine[i] = rasterParams[1];
+  }
+  else
+  {
+    // 180.0 / M_PI = 57.29577951308232
+    float aspect = atan2( derX, derY ) * 57.29577951308232;
+    if ( aspect < 0 )
+        resultLine[i] = 90.0f - aspect;
+    else if (aspect > 90.0f)
+        // 360 + 90 = 450
+        resultLine[i] = 450.0f - aspect;
+    else
+        resultLine[i] = 90.0 - aspect;
+  }
+}

resources/opencl_programs/calcfirstder.cl

@@ -0,0 +1,71 @@
+// Calculate the first derivative from a 3x3 cell matrix
+float calcFirstDer( float x11, float x21, float x31, float x12, float x22, float x32, float x13, float x23, float x33,
+                       float inputNodataValue, float outputNodataValue, float zFactor, float mCellSize )
+{
+ //the basic formula would be simple, but we need to test for nodata values...
+ //X: return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * cellSizeX));
+ //Y: return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * cellSizeY));
+
+ int weight = 0;
+ float sum = 0;
+
+
+ //first row
+ if ( x31 != inputNodataValue && x11 != inputNodataValue ) //the normal case
+ {
+   sum += ( x31 - x11 );
+   weight += 2;
+ }
+ else if ( x31 == inputNodataValue && x11 != inputNodataValue && x21 != inputNodataValue ) //probably 3x3 window is at the border
+ {
+   sum += ( x21 - x11 );
+   weight += 1;
+ }
+ else if ( x11 == inputNodataValue && x31 != inputNodataValue && x21 != inputNodataValue ) //probably 3x3 window is at the border
+ {
+   sum += ( x31 - x21 );
+   weight += 1;
+ }
+
+ //second row
+ if ( x32 != inputNodataValue && x12 != inputNodataValue ) //the normal case
+ {
+   sum += 2.0f * ( x32 - x12 );
+   weight += 4;
+ }
+ else if ( x32 == inputNodataValue && x12 != inputNodataValue && x22 != inputNodataValue )
+ {
+   sum += 2.0f * ( x22 - x12 );
+   weight += 2;
+ }
+ else if ( x12 == inputNodataValue && x32 != inputNodataValue && x22 != inputNodataValue )
+ {
+   sum += 2.0f * ( x32 - x22 );
+   weight += 2;
+ }
+
+ //third row
+ if ( x33 != inputNodataValue && x13 != inputNodataValue ) //the normal case
+ {
+   sum += ( x33 - x13 );
+   weight += 2;
+ }
+ else if ( x33 == inputNodataValue && x13 != inputNodataValue && x23 != inputNodataValue )
+ {
+   sum += ( x23 - x13 );
+   weight += 1;
+ }
+ else if ( x13 == inputNodataValue && x33 != inputNodataValue && x23 != inputNodataValue )
+ {
+   sum += ( x33 - x23 );
+   weight += 1;
+ }
+
+ if ( weight == 0 )
+ {
+   return outputNodataValue;
+ }
+
+ return sum / ( weight * mCellSize ) * zFactor;
+}
+

resources/opencl_programs/slope.cl

@@ -0,0 +1,38 @@
+#include "calcfirstder.cl"
+
+__kernel void processNineCellWindow( __global float *scanLine1,
+                                     __global float *scanLine2,
+                                     __global float *scanLine3,
+                                     __global float *resultLine,
+                                     __global float *rasterParams
+                       ) {
+
+  // Get the index of the current element
+  const int i = get_global_id(0);
+
+  // Do the operation
+  //return (( (x31 - x11) + 2 * (x32 - x12) + (x33 - x13) ) / (8 * mCellSizeX))
+  float derX = calcFirstDer(   scanLine1[i],   scanLine2[i],   scanLine3[i],
+                               scanLine1[i+1], scanLine2[i+1], scanLine3[i+1],
+                               scanLine1[i+2], scanLine2[i+2], scanLine3[i+2],
+                               rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3]
+                             );
+  //return (((x11 - x13) + 2 * (x21 - x23) + (x31 - x33)) / ( 8 * mCellSizeY));
+  float derY = calcFirstDer(   scanLine1[i+2], scanLine1[i+1], scanLine1[i],
+                               scanLine2[i+2], scanLine2[i+1], scanLine2[i],
+                               scanLine3[i+2], scanLine3[i+1], scanLine3[i],
+                               rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]
+                             );
+
+
+  if ( derX == rasterParams[1] || derY == rasterParams[1] )
+  {
+    resultLine[i] = rasterParams[1];
+  }
+  else
+  {
+    float res = sqrt( derX * derX + derY * derY );
+    res = atanpi( res );
+    resultLine[i] = res * 180.0;
+  }
+}

src/analysis/raster/qgsninecellfilter.cpp

@@ -41,7 +41,7 @@ int QgsNineCellFilter::processRaster( QgsFeedback *feedback )
   if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! openClProgramBaseName( ).isEmpty() )
   {
     // Load the program sources
-    QString source( QgsOpenClUtils::sourceFromPath( QStringLiteral( "/home/ale/dev/QGIS/src/analysis/raster/%1.cl" ).arg( openClProgramBaseName( ) ) ) );
+    QString source( QgsOpenClUtils::sourceFromBaseName( openClProgramBaseName( ) ) );
     if ( ! source.isEmpty() )
     {
       try
@@ -50,19 +50,6 @@ int QgsNineCellFilter::processRaster( QgsFeedback *feedback )
                                    .arg( openClProgramBaseName( ) ), QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Info );
         return processRasterGPU( source, feedback );
       }
-      catch ( cl::BuildError &e )
-      {
-        cl::BuildLogType build_logs = e.getBuildLog();
-        QString build_log;
-        if ( build_logs.size() > 0 )
-          build_log = QString::fromStdString( build_logs[0].second );
-        else
-          build_log = QObject::tr( "Build logs not available!" );
-        QString err = QObject::tr( "Error building OpenCL program: %1" )
-                      .arg( build_log );
-        QgsMessageLog::logMessage( err, QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
-        throw QgsProcessingException( err );
-      }
       catch ( cl::Error &e )
       {
         QString err = QObject::tr( "Error %1 running OpenCL program in %2" )
@@ -220,18 +207,15 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
     return 6;
   }
 
-  // Prepare context
+  // Prepare context and queue
   cl::Context ctx = QgsOpenClUtils::context();
   cl::Context::setDefault( ctx );
+  cl::CommandQueue queue( ctx );
 
   //keep only three scanlines in memory at a time, make room for initial and final nodata
   QgsOpenClUtils::CPLAllocator<float> scanLine1( xSize + 2 );
   QgsOpenClUtils::CPLAllocator<float> scanLine2( xSize + 2 );
   QgsOpenClUtils::CPLAllocator<float> scanLine3( xSize + 2 );
-  //float *scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * ( xSize + 2 ) );
-  //float *scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * ( xSize + 2 ) );
-
-  //float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
   QgsOpenClUtils::CPLAllocator<float> resultLine( xSize );
 
   cl_int errorCode = 0;
@@ -255,9 +239,7 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
   cl::Buffer resultLineBuffer( CL_MEM_WRITE_ONLY, sizeof( float ) * xSize, nullptr, &errorCode );
 
   // Create a program from the kernel source
-  cl::Program program( source.toStdString() );
-  // Use CL 1.1 for compatibility with older libs
-  program.build( "-cl-std=CL1.1" );
+  cl::Program program( QgsOpenClUtils::buildProgram( ctx, source, QgsOpenClUtils::ExceptionBehavior::Throw ) );
 
   // Create the OpenCL kernel
   auto kernel = cl::KernelFunctor <
@@ -297,9 +279,6 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
     else
     {
       //normally fetch only scanLine3 and release scanline 1 if we move forward one row
-      //scanLine1 = scanLine2;
-      //scanLine2 = scanLine3;
-      //scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * ( xSize + 2 ) );
       scanLine1.reset( scanLine2.release() );
       scanLine2.reset( scanLine3.release() );
       scanLine3.reset( xSize + 2 );
@@ -325,6 +304,9 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
     scanLine2[0] = scanLine2[xSize + 1] = mInputNodataValue;
     scanLine3[0] = scanLine3[xSize + 1] = mInputNodataValue;
 
+    // TODO: There is room for further optimization here: instead of replacing the buffers
+    //       we could just replace just hthe new one (the top row) and switch the order
+    //       of buffer arguments in the kernell call.
     errorCode = cl::enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0,
                                         sizeof( float ) * ( xSize + 2 ), scanLine1.get() );
     errorCode = cl::enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0,
@@ -333,6 +315,7 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
                                         sizeof( float ) * ( xSize + 2 ), scanLine3.get() );
 
     kernel( cl::EnqueueArgs(
+              queue,
               cl::NDRange( xSize )
             ),
             scanLine1Buffer,