Merge pull request #7451 from elpaso/opencl-utils-2

[feature] OpenCL support
qgis · Aug 8, 2018 · 0b502ff · 0b502ff
2 parents 55473e7 + 0960b1f
commit 0b502ff
Show file tree

Hide file tree

Showing 45 changed files with 2,752 additions and 331 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -28,6 +28,18 @@ MESSAGE(STATUS "QGIS version: ${COMPLETE_VERSION} ${RELEASE_NAME} (${QGIS_VERSIO
 
 
 #############################################################
+# Configure OpenCL if available
+
+OPTION(USE_OPENCL "Use OpenCL" ON)
+IF (USE_OPENCL)
+    FIND_PACKAGE(OpenCL)
+    IF(${OpenCL_FOUND})
+        SET(HAVE_OPENCL TRUE)
+    ELSE(${OpenCL_FOUND})
+        MESSAGE(STATUS "Couldn't find OpenCL: support DISABLED")
+    ENDIF(${OpenCL_FOUND})
+ENDIF(USE_OPENCL)
+
 # Configure CCache if available
 IF(NOT MSVC)
   option(USE_CCACHE "Use ccache" ON)

diff --git a/cmake_templates/qgsconfig.h.in b/cmake_templates/qgsconfig.h.in
@@ -57,6 +57,8 @@
 
 #cmakedefine HAVE_SERVER_PYTHON_PLUGINS
 
+#cmakedefine HAVE_OPENCL
+
 #cmakedefine ENABLE_MODELTEST
 
 #cmakedefine HAVE_3D

diff --git a/images/images.qrc b/images/images.qrc
@@ -699,6 +699,7 @@
         <file>themes/default/mIndicatorEmbedded.svg</file>
         <file>themes/default/mIconHistory.svg</file>
         <file>themes/default/mIndicatorMemory.svg</file>
+        <file>themes/default/mIconGPU.svg</file>
     </qresource>
     <qresource prefix="/images/tips">
         <file alias="symbol_levels.png">qgis_tips/symbol_levels.png</file>

diff --git a/images/themes/default/mIconGPU.svg b/images/themes/default/mIconGPU.svg
diff --git a/python/analysis/auto_generated/raster/qgsaspectfilter.sip.in b/python/analysis/auto_generated/raster/qgsaspectfilter.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
 
+
 };
 
 /************************************************************************

diff --git a/python/analysis/auto_generated/raster/qgsninecellfilter.sip.in b/python/analysis/auto_generated/raster/qgsninecellfilter.sip.in
@@ -10,6 +10,7 @@
 
 
 
+
 class QgsNineCellFilter
 {
 %Docstring
@@ -34,7 +35,7 @@ Starts the calculation, reads from mInputFile and stores the result in mOutputFi
 
 :param feedback: feedback object that receives update and that is checked for cancelation.
 
-:return: 0 in case of success*
+:return: 0 in case of success
 %End
 
     double cellSizeX() const;
@@ -54,8 +55,23 @@ Starts the calculation, reads from mInputFile and stores the result in mOutputFi
                                          float *x12, float *x22, float *x32,
                                          float *x13, float *x23, float *x33 ) = 0;
 %Docstring
-Calculates output value from nine input values. The input values and the output value can be equal to the
-nodata value if not present or outside of the border. Must be implemented by subclasses*
+Calculates output value from nine input values. The input values and the output
+value can be equal to the nodata value if not present or outside of the border.
+Must be implemented by subclasses.
+
+First index of the input cell is the row, second index is the column
+
+:param x11: surrounding cell top left
+:param x21: surrounding cell central left
+:param x31: surrounding cell bottom left
+:param x12: surrounding cell top central
+:param x22: the central cell for which the value will be calculated
+:param x32: surrounding cell bottom central
+:param x13: surrounding cell top right
+:param x23: surrounding cell central right
+:param x33: surrounding cell bottom right
+
+:return: the calculated cell value for the central cell x22
 %End
 
   protected:

diff --git a/python/analysis/auto_generated/raster/qgsslopefilter.sip.in b/python/analysis/auto_generated/raster/qgsslopefilter.sip.in
@@ -28,6 +28,8 @@ Calculates slope values in a window of 3x3 cells based on first order derivative
 Calculates output value from nine input values. The input values and the output value can be equal to the
 nodata value if not present or outside of the border. Must be implemented by subclasses*
 %End
+
+
 };
 
 /************************************************************************

diff --git a/python/core/auto_generated/raster/qgsrasterblock.sip.in b/python/core/auto_generated/raster/qgsrasterblock.sip.in
@@ -186,7 +186,7 @@ Read a single value
 :return: color *
 %End
 
-    bool isNoData( int row, int column );
+    bool isNoData( int row, int column ) const;
 %Docstring
 Check if value at position is no data
 
@@ -196,7 +196,17 @@ Check if value at position is no data
 :return: true if value is no data *
 %End
 
-    bool isNoData( qgssize index );
+    bool isNoData( qgssize row, qgssize column ) const;
+%Docstring
+Check if value at position is no data
+
+:param row: row index
+:param column: column index
+
+:return: true if value is no data *
+%End
+
+    bool isNoData( qgssize index ) const;
 %Docstring
 Check if value at position is no data
 

diff --git a/resources/CMakeLists.txt b/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)

diff --git a/resources/opencl_programs/aspect.cl b/resources/opencl_programs/aspect.cl
@@ -0,0 +1,40 @@
+#include "calcfirstder.cl"
+
+__kernel void processNineCellWindow( __global float *scanLine1,
+                                     __global float *scanLine2,
+                                     __global float *scanLine3,
+                                     __global float *resultLine,
+                                     __global float *rasterParams // mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY
+                                     )
+{
+
+  // Get the index of the current element
+  const int i = get_global_id(0);
+
+  if ( scanLine2[i+1] == rasterParams[0] )
+  {
+     resultLine[i] = rasterParams[1];
+  }
+  else
+  {
+    float derX = calcFirstDer( scanLine1[i], scanLine1[i+1], scanLine1[i+2],
+                               scanLine2[i], scanLine2[i+1], scanLine2[i+2],
+                               scanLine3[i], scanLine3[i+1], scanLine3[i+2],
+                               rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3] );
+
+    float derY = calcFirstDer( scanLine3[i],   scanLine2[i],   scanLine1[i],
+                               scanLine3[i+1], scanLine2[i+1], scanLine1[i+1],
+                               scanLine3[i+2], scanLine2[i+2], scanLine1[i+2],
+                               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
+    {
+      resultLine[i] = 180.0f + atan2pi( derX, derY ) * 180.0f;
+    }
+  }
+}
diff --git a/resources/opencl_programs/aspect_renderer.cl b/resources/opencl_programs/aspect_renderer.cl
@@ -0,0 +1,53 @@
+#include "calcfirstder.cl"
+
+// Aspect renderer for QGIS
+
+__kernel void processNineCellWindow( __global float *scanLine1,
+                                     __global float *scanLine2,
+                                     __global float *scanLine3,
+                                     __global uchar4 *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]
+                             );
+
+
+  float res;
+  if ( derX == rasterParams[1] || derY == rasterParams[1] ||
+       ( derX == 0.0f && derY == 0.0f) )
+  {
+    res = rasterParams[1];
+  }
+  else
+  {
+    // 180.0 / M_PI = 57.29577951308232
+    float aspect = atan2( derX, derY ) * 57.29577951308232;
+    if ( aspect < 0 )
+        res = 90.0f - aspect;
+    else if (aspect > 90.0f)
+        // 360 + 90 = 450
+        res = 450.0f - aspect;
+    else
+        res = 90.0 - aspect;
+  }
+
+  res = res / 360 * 255;
+
+  resultLine[i] = (uchar4)(res, res, res, 255);
+}
+
diff --git a/resources/opencl_programs/calcfirstder.cl b/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;
+}
+
diff --git a/resources/opencl_programs/hillshade.cl b/resources/opencl_programs/hillshade.cl
@@ -0,0 +1,56 @@
+#include "calcfirstder.cl"
+
+#ifndef M_PI
+    #define M_PI 3.14159265358979323846
+#endif
+
+__kernel void processNineCellWindow( __global float *scanLine1,
+                                     __global float *scanLine2,
+                                     __global float *scanLine3,
+                                     __global float *resultLine,
+                                     __global float *rasterParams // mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY, zenith_rad, azimuth_rad
+                                     )
+{
+
+    // Get the index of the current element
+    const int i = get_global_id(0);
+
+    if ( scanLine2[i+1] == rasterParams[0] )
+    {
+       resultLine[i] = rasterParams[1];
+    }
+    else
+    {
+
+      float derX = calcFirstDer( scanLine1[i], scanLine1[i+1], scanLine1[i+2],
+                                 scanLine2[i], scanLine2[i+1], scanLine2[i+2],
+                                 scanLine3[i], scanLine3[i+1], scanLine3[i+2],
+                                 rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[3] );
+
+      float derY = calcFirstDer( scanLine3[i],   scanLine2[i],   scanLine1[i],
+                                 scanLine3[i+1], scanLine2[i+1], scanLine1[i+1],
+                                 scanLine3[i+2], scanLine2[i+2], scanLine1[i+2],
+                                 rasterParams[0], rasterParams[1], rasterParams[2], rasterParams[4]);
+
+      if ( derX == rasterParams[1] || derY == rasterParams[1] )
+      {
+        resultLine[i] = rasterParams[1];
+      }
+      else
+      {
+
+        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 = rasterParams[7] / 2.0f;
+        else
+          aspect_rad = M_PI + atan2( derX, derY );
+
+        resultLine[i] = max(0.0f, 255.0f * ( ( rasterParams[5] * cos( slope_rad ) ) +
+                                            (  rasterParams[6] * sin( slope_rad ) *
+                                              cos( rasterParams[7] - aspect_rad ) ) ) );
+
+      }
+    }
+}