diff --git a/CMakeLists.txt b/CMakeLists.txt index 22bbcae9f7f..074afdca60b 100644 --- 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 index 2310ac83d9e..e267d32e5af 100644 --- 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 index a622a00af06..cbb09505033 100755 --- a/images/images.qrc +++ b/images/images.qrc @@ -699,6 +699,7 @@ themes/default/mIndicatorEmbedded.svg themes/default/mIconHistory.svg themes/default/mIndicatorMemory.svg + themes/default/mIconGPU.svg qgis_tips/symbol_levels.png diff --git a/images/themes/default/mIconGPU.svg b/images/themes/default/mIconGPU.svg new file mode 100644 index 00000000000..898b2f5dcc9 --- /dev/null +++ b/images/themes/default/mIconGPU.svg @@ -0,0 +1,14 @@ + + + + + + diff --git a/python/analysis/auto_generated/raster/qgsaspectfilter.sip.in b/python/analysis/auto_generated/raster/qgsaspectfilter.sip.in index fb5829bd81b..d9b89d526b2 100644 --- 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 index 023cc8625b9..87eb133f2ee 100644 --- 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 index 3f5f9a7e244..59a0de0ee5b 100644 --- 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 index b81f4991730..cc2e6e79b88 100644 --- 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 index 32e09478ee6..f1e755c5987 100644 --- 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 new file mode 100644 index 00000000000..91b829e0d84 --- /dev/null +++ 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 new file mode 100644 index 00000000000..7de6401dbe3 --- /dev/null +++ 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 new file mode 100644 index 00000000000..e23caf123a4 --- /dev/null +++ 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 new file mode 100644 index 00000000000..000d0700a07 --- /dev/null +++ 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 ) ) ) ); + + } + } +} diff --git a/resources/opencl_programs/hillshade_renderer.cl b/resources/opencl_programs/hillshade_renderer.cl new file mode 100644 index 00000000000..e0d25bc63e6 --- /dev/null +++ b/resources/opencl_programs/hillshade_renderer.cl @@ -0,0 +1,110 @@ +__kernel void processNineCellWindow( __global float *scanLine1, + __global float *scanLine2, + __global float *scanLine3, + __global uchar4 *resultLine, // This is an image! + __global float *rasterParams + + ) +{ + + // Get the index of the current element + const int i = get_global_id(0); + + float x11 = scanLine1[i]; + float x12 = scanLine1[i+1]; + float x13 = scanLine1[i+2]; + float x21 = scanLine2[i]; + float x22 = scanLine2[i+1]; + float x23 = scanLine2[i+2]; + float x31 = scanLine3[i]; + float x32 = scanLine3[i+1]; + float x33 = scanLine3[i+2]; + + if ( x22 == rasterParams[0] ) + { + float alpha = rasterParams[8] * rasterParams[16]; + resultLine[i] = (uchar4)(rasterParams[13] * alpha, + rasterParams[14] * alpha, + rasterParams[15] * alpha, 255 * alpha); + } + else + { + if ( x11 == rasterParams[0] ) x11 = x22; + if ( x12 == rasterParams[0] ) x12 = x22; + if ( x13 == rasterParams[0] ) x13 = x22; + if ( x21 == rasterParams[0] ) x21 = x22; + // Note: skip central cell x22 + if ( x23 == rasterParams[0] ) x23 = x22; + if ( x31 == rasterParams[0] ) x31 = x22; + if ( x32 == rasterParams[0] ) x32 = x22; + if ( x33 == rasterParams[0] ) x33 = x22; + + float derX = ( ( x13 + x23 + x23 + x33 ) - ( x11 + x21 + x21 + x31 ) ) / ( 8 * rasterParams[3] ); + float derY = ( ( x31 + x32 + x32 + x33 ) - ( x11 + x12 + x12 + x13 ) ) / ( 8 * -rasterParams[4]); + + if ( derX == rasterParams[0] || + derX == rasterParams[0] ) + { + float alpha = rasterParams[5] * rasterParams[16]; + resultLine[i] = (uchar4)(rasterParams[13] * alpha, + rasterParams[14] * alpha, + rasterParams[15] * alpha, 255 * alpha); + } + else + { + float res; + if ( rasterParams[17] ) // Multi directional? + { + // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf + // Fast formula from GDAL DEM + const float xx = derX * derX; + const float yy = derY * derY; + const float xx_plus_yy = xx + yy; + // Flat? + if ( xx_plus_yy == 0.0f ) + { + res = clamp( 1.0f + rasterParams[9], 0.0f, 255.0f ); + } + else + { + // ... then the shade value from different azimuth + float val225_mul_127 = rasterParams[10] + + ( derX - derY ) * rasterParams[11]; + val225_mul_127 = ( val225_mul_127 <= 0.0f ) ? 0.0f : val225_mul_127; + float val270_mul_127 = rasterParams[10] - + derX * rasterParams[12]; + val270_mul_127 = ( val270_mul_127 <= 0.0f ) ? 0.0f : val270_mul_127; + float val315_mul_127 = rasterParams[10] + + ( derX + derY ) * rasterParams[11]; + val315_mul_127 = ( val315_mul_127 <= 0.0f ) ? 0.0f : val315_mul_127; + float val360_mul_127 = rasterParams[10] - + derY * rasterParams[12]; + val360_mul_127 = ( val360_mul_127 <= 0.0f ) ? 0.0f : val360_mul_127; + // ... then the weighted shading + const float weight_225 = 0.5f * xx_plus_yy - derX * derY; + const float weight_270 = xx; + const float weight_315 = xx_plus_yy - weight_225; + const float weight_360 = yy; + const float cang_mul_127 = ( + ( weight_225 * val225_mul_127 + + weight_270 * val270_mul_127 + + weight_315 * val315_mul_127 + + weight_360 * val360_mul_127 ) / xx_plus_yy ) / + ( 1 + rasterParams[8] * xx_plus_yy ); + res = clamp( 1.0f + cang_mul_127, 0.0f, 255.0f ); + } + } + else + { + res = ( rasterParams[9] - + ( derY * rasterParams[6] - + derX * rasterParams[7] )) / + sqrt( 1 + rasterParams[8] * + ( derX * derX + derY * derY ) ); + res = res <= 0.0f ? 1.0f : 1.0f + res; + } + res = res * rasterParams[5]; + resultLine[i] = (uchar4)(res, res, res, 255 * rasterParams[5]); + } + } +} diff --git a/resources/opencl_programs/ruggedness.cl b/resources/opencl_programs/ruggedness.cl new file mode 100644 index 00000000000..396a3467b47 --- /dev/null +++ b/resources/opencl_programs/ruggedness.cl @@ -0,0 +1,51 @@ +__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); + + float x11 = scanLine1[i]; + float x21 = scanLine1[i+1]; + float x31 = scanLine1[i+2]; + float x12 = scanLine2[i]; + float x22 = scanLine2[i+1]; + float x32 = scanLine2[i+2]; + float x13 = scanLine3[i]; + float x23 = scanLine3[i+1]; + float x33 = scanLine3[i+2]; + + if ( x22 == rasterParams[0] ) + { + resultLine[i] = rasterParams[1]; + } + else + { + // Nodata handling + if ( x11 == rasterParams[0] ) x11 = x22; + if ( x12 == rasterParams[0] ) x12 = x22; + if ( x13 == rasterParams[0] ) x13 = x22; + if ( x21 == rasterParams[0] ) x21 = x22; + if ( x23 == rasterParams[0] ) x23 = x22; + if ( x31 == rasterParams[0] ) x31 = x22; + if ( x32 == rasterParams[0] ) x32 = x22; + if ( x33 == rasterParams[0] ) x33 = x22; + + float diff1 = x11 - x22; + float diff2 = x21 - x22; + float diff3 = x31 - x22; + float diff4 = x12 - x22; + float diff5 = x32 - x22; + float diff6 = x13 - x22; + float diff7 = x23 - x22; + float diff8 = x33 - x22; + + resultLine[i] = sqrt(diff1 * diff1 + diff2 * diff2 + + diff3 * diff3 + diff4 * diff4 + + diff5 * diff5 + diff6 * diff6 + + diff7 * diff7 + diff8 * diff8); + } +} diff --git a/resources/opencl_programs/slope.cl b/resources/opencl_programs/slope.cl new file mode 100644 index 00000000000..3113970ad63 --- /dev/null +++ b/resources/opencl_programs/slope.cl @@ -0,0 +1,42 @@ +#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] ) + { + resultLine[i] = rasterParams[1]; + } + else + { + float res = sqrt( derX * derX + derY * derY ); + res = atanpi( res ); + resultLine[i] = res * 180.0f; + } + } +} diff --git a/resources/opencl_programs/slope_renderer.cl b/resources/opencl_programs/slope_renderer.cl new file mode 100644 index 00000000000..674c7ca6a02 --- /dev/null +++ b/resources/opencl_programs/slope_renderer.cl @@ -0,0 +1,45 @@ +#include "calcfirstder.cl" + +// Slope renderer for QGIS + +__kernel void processNineCellWindow( __global float *scanLine1, + __global float *scanLine2, + __global float *scanLine3, + __global uchar4 *resultLine, // This is an image BGRA ! + __global float *rasterParams // [ mInputNodataValue, mOutputNodataValue, mZFactor, mCellSizeX, mCellSizeY ] + + ) { + + // 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] ) + { + res = rasterParams[1]; + } + else + { + float slope = sqrt( derX * derX + derY * derY ); + slope = atan( slope ); + res = slope * 255; + } + + resultLine[i] = (uchar4)(res, res, res, 255); + +} diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index 244cbf90fad..ac188e47067 100755 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -302,6 +302,12 @@ INCLUDE_DIRECTORIES(SYSTEM ${SQLITE3_INCLUDE_DIR} ) +IF(HAVE_OPENCL) + INCLUDE_DIRECTORIES(${OpenCL_INCLUDE_DIRS}) +ENDIF(HAVE_OPENCL) + + + ############################################################# # qgis_analysis library @@ -353,6 +359,11 @@ ENDIF (NOT ANDROID) TARGET_LINK_LIBRARIES(qgis_analysis qgis_core) +IF(HAVE_OPENCL) + TARGET_LINK_LIBRARIES(qgis_analysis ${OpenCL_LIBRARIES}) +ENDIF(HAVE_OPENCL) + + # clang-tidy IF(CLANG_TIDY_EXE) SET_TARGET_PROPERTIES( diff --git a/src/analysis/raster/qgsalignraster.cpp b/src/analysis/raster/qgsalignraster.cpp index cacc90b6b0d..60cd52366d0 100644 --- a/src/analysis/raster/qgsalignraster.cpp +++ b/src/analysis/raster/qgsalignraster.cpp @@ -615,3 +615,4 @@ double QgsAlignRaster::RasterInfo::identify( double mx, double my ) return value; } + diff --git a/src/analysis/raster/qgsaspectfilter.h b/src/analysis/raster/qgsaspectfilter.h index d04ecf442a7..8d271ad85bd 100644 --- a/src/analysis/raster/qgsaspectfilter.h +++ b/src/analysis/raster/qgsaspectfilter.h @@ -36,6 +36,16 @@ class ANALYSIS_EXPORT QgsAspectFilter: public QgsDerivativeFilter float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) override; + +#ifdef HAVE_OPENCL + private: + + const QString openClProgramBaseName() const override + { + return QStringLiteral( "aspect" ); + } +#endif + }; #endif // QGSASPECTFILTER_H diff --git a/src/analysis/raster/qgshillshadefilter.cpp b/src/analysis/raster/qgshillshadefilter.cpp index b0ee79da57c..90cafa37dcb 100644 --- a/src/analysis/raster/qgshillshadefilter.cpp +++ b/src/analysis/raster/qgshillshadefilter.cpp @@ -21,8 +21,11 @@ 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( lightAzimuth ) ) + , mLightAngle( static_cast( lightAngle ) ) + , mCosZenithRad( std::cos( static_cast( lightAngle * M_PI ) / 180.0f ) ) + , mSinZenithRad( std::sin( static_cast( lightAngle * M_PI ) / 180.0f ) ) + , mAzimuthRad( static_cast( lightAzimuth * M_PI ) / 180.0f ) { } @@ -30,6 +33,7 @@ float QgsHillshadeFilter::processNineCellWindow( float *x11, float *x21, float * 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( M_PI ) / 180.0f; +} + +void QgsHillshadeFilter::setLightAngle( float angle ) +{ + mLightAngle = angle; + mCosZenithRad = std::cos( angle * static_cast( M_PI ) / 180.0f ); + mSinZenithRad = std::sin( angle * static_cast( M_PI ) / 180.0f ); +} + +#ifdef HAVE_OPENCL + +void QgsHillshadeFilter::addExtraRasterParams( std::vector ¶ms ) +{ + + params.push_back( mCosZenithRad ); // cos_zenith_rad 5 + params.push_back( mSinZenithRad ); // sin_zenith_rad 6 + params.push_back( mAzimuthRad ); // azimuth_rad 7 + +} + +#endif diff --git a/src/analysis/raster/qgshillshadefilter.h b/src/analysis/raster/qgshillshadefilter.h index 8a5069d9faf..86d03d63583 100644 --- a/src/analysis/raster/qgshillshadefilter.h +++ b/src/analysis/raster/qgshillshadefilter.h @@ -39,13 +39,34 @@ 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 ¶ms ) override; +#endif + }; #endif // QGSHILLSHADEFILTER_H diff --git a/src/analysis/raster/qgsninecellfilter.cpp b/src/analysis/raster/qgsninecellfilter.cpp index ae0054d705a..4c185775120 100644 --- a/src/analysis/raster/qgsninecellfilter.cpp +++ b/src/analysis/raster/qgsninecellfilter.cpp @@ -20,7 +20,18 @@ #include "cpl_string.h" #include "qgsfeedback.h" #include "qgsogrutils.h" +#include "qgsmessagelog.h" + +#ifdef HAVE_OPENCL +#include "qgsopenclutils.h" +#endif + #include +#include +#include +#include + + QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat ) : mInputFile( inputFile ) @@ -30,146 +41,44 @@ QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &o } +// TODO: return an anum instead of an int int QgsNineCellFilter::processRaster( QgsFeedback *feedback ) { - GDALAllRegister(); - - //open input file - int xSize, ySize; - gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) ); - if ( !inputDataset ) +#ifdef HAVE_OPENCL + if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! openClProgramBaseName( ).isEmpty() ) { - return 1; //opening of input file failed - } - - //output driver - GDALDriverH outputDriver = openOutputDriver(); - if ( !outputDriver ) - { - return 2; - } - - gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) ); - if ( !outputDataset ) - { - return 3; //create operation on output file failed - } - - //open first raster band for reading (operation is only for single band raster) - GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 ); - if ( !rasterBand ) - { - return 4; - } - mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr ); - - GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 ); - if ( !outputRasterBand ) - { - return 5; - } - //try to set -9999 as nodata value - GDALSetRasterNoDataValue( outputRasterBand, -9999 ); - mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr ); - - if ( ySize < 3 ) //we require at least three rows (should be true for most datasets) - { - return 6; - } - - //keep only three scanlines in memory at a time - float *scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); - float *scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); - float *scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); - - 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 ) - { - if ( feedback && feedback->isCanceled() ) + // Load the program sources + QString source( QgsOpenClUtils::sourceFromBaseName( openClProgramBaseName( ) ) ); + if ( ! source.isEmpty() ) { - break; - } - - if ( feedback ) - { - feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize ); - } - - if ( i == 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; ++a ) + try { - scanLine1[a] = mInputNodataValue; + QgsDebugMsg( QObject::tr( "Running OpenCL program: %1" ).arg( openClProgramBaseName( ) ) ); + return processRasterGPU( source, feedback ); } - if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + catch ( cl::Error &e ) { - QgsDebugMsg( "Raster IO Error" ); + QString err = QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what( ) ).arg( QgsOpenClUtils::errorText( e.err( ) ) ); + QgsMessageLog::logMessage( err, QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical ); + throw QgsProcessingException( err ); } } else { - //normally fetch only scanLine3 and release scanline 1 if we move forward one row - CPLFree( scanLine1 ); - scanLine1 = scanLine2; - scanLine2 = scanLine3; - scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); - } - - if ( i == ySize - 1 ) //fill the row below the bottom with nodata values - { - for ( int a = 0; a < xSize; ++a ) - { - scanLine3[a] = mInputNodataValue; - } - } - else - { - if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) - { - QgsDebugMsg( "Raster IO Error" ); - } - } - - for ( int j = 0; j < xSize; ++j ) - { - if ( j == 0 ) - { - resultLine[j] = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j + 1], &mInputNodataValue, &scanLine2[j], - &scanLine2[j + 1], &mInputNodataValue, &scanLine3[j], &scanLine3[j + 1] ); - } - else if ( j == xSize - 1 ) - { - resultLine[j] = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &mInputNodataValue, &scanLine2[j - 1], &scanLine2[j], - &mInputNodataValue, &scanLine3[j - 1], &scanLine3[j], &mInputNodataValue ); - } - else - { - resultLine[j] = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &scanLine1[j + 1], &scanLine2[j - 1], &scanLine2[j], - &scanLine2[j + 1], &scanLine3[j - 1], &scanLine3[j], &scanLine3[j + 1] ); - } - } - - if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) - { - QgsDebugMsg( "Raster IO Error" ); + QString err = QObject::tr( "Error loading OpenCL program sources" ); + QgsMessageLog::logMessage( err, + QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical ); + throw QgsProcessingException( err ); } } - - CPLFree( resultLine ); - CPLFree( scanLine1 ); - CPLFree( scanLine2 ); - CPLFree( scanLine3 ); - - if ( feedback && feedback->isCanceled() ) + else { - //delete the dataset without closing (because it is faster) - gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile ); - return 7; + return processRasterCPU( feedback ); } - return 0; + return 1; +#else + return processRasterCPU( feedback ); +#endif } gdal::dataset_unique_ptr QgsNineCellFilter::openInputFile( int &nCellsX, int &nCellsY ) @@ -254,3 +163,334 @@ gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDa return outputDataset; } +#ifdef HAVE_OPENCL + +// TODO: return an anum instead of an int +int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *feedback ) +{ + + GDALAllRegister(); + + //open input file + int xSize, ySize; + gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) ); + if ( !inputDataset ) + { + return 1; //opening of input file failed + } + + //output driver + GDALDriverH outputDriver = openOutputDriver(); + if ( !outputDriver ) + { + return 2; + } + + gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) ); + if ( !outputDataset ) + { + return 3; //create operation on output file failed + } + + //open first raster band for reading (operation is only for single band raster) + GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 ); + if ( !rasterBand ) + { + return 4; + } + mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr ); + + GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 ); + if ( !outputRasterBand ) + { + return 5; + } + //try to set -9999 as nodata value + GDALSetRasterNoDataValue( outputRasterBand, -9999 ); + mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr ); + + if ( ySize < 3 ) //we require at least three rows (should be true for most datasets) + { + return 6; + } + + // Prepare context and queue + cl::Context ctx = QgsOpenClUtils::context(); + cl::CommandQueue queue( ctx ); + + //keep only three scanlines in memory at a time, make room for initial and final nodata + QgsOpenClUtils::CPLAllocator scanLine( xSize + 2 ); + QgsOpenClUtils::CPLAllocator resultLine( xSize ); + + // Cast to float (because double just crashes on some GPUs) + std::vector rasterParams; + + rasterParams.push_back( mInputNodataValue ); // 0 + rasterParams.push_back( mOutputNodataValue ); // 1 + rasterParams.push_back( mZFactor ); // 2 + rasterParams.push_back( mCellSizeX ); // 3 + rasterParams.push_back( mCellSizeY ); // 4 + + // Allow subclasses to add extra params needed for computation: + // used to pass additional args to opencl program + addExtraRasterParams( rasterParams ); + + std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) ); + std::size_t inputSize( sizeof( float ) * ( xSize ) ); + + cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr ); + cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr ); + cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr ); + cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr ); + cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer}; + cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, inputSize, nullptr, nullptr ); + + // Create a program from the kernel source + cl::Program program( QgsOpenClUtils::buildProgram( ctx, source, QgsOpenClUtils::ExceptionBehavior::Throw ) ); + + // Create the OpenCL kernel + auto kernel = cl::KernelFunctor < + cl::Buffer &, + cl::Buffer &, + cl::Buffer &, + cl::Buffer &, + cl::Buffer & + > ( program, "processNineCellWindow" ); + + // Rotate buffer index + std::vector rowIndex = {0, 1, 2}; + + // 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 ) + { + if ( feedback && feedback->isCanceled() ) + { + break; + } + + if ( feedback ) + { + feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize ); + } + + if ( i == 0 ) + { + // Fill scanline 1 with (input) nodata for the values above the first row and + // feed scanline2 with the first actual data row + for ( int a = 0; a < xSize + 2 ; ++a ) + { + scanLine[a] = mInputNodataValue; + } + queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine.get() ); + + // Read scanline2: first real raster row + if ( GDALRasterIO( rasterBand, GF_Read, 0, i, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + QgsDebugMsg( "Raster IO Error" ); + } + queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine.get() ); + + // Read scanline3: second real raster row + if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + QgsDebugMsg( "Raster IO Error" ); + } + queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine.get() ); + } + else + { + // Normally fetch only scanLine3 and move forward one row + // Read scanline 3, fill the last row with nodata values if it's the last iteration + if ( i == ySize - 1 ) //fill the row below the bottom with nodata values + { + for ( int a = 0; a < xSize + 2; ++a ) + { + scanLine[a] = mInputNodataValue; + } + queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0 + } + 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 ) + { + QgsDebugMsg( "Raster IO Error" ); + } + queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0 + } + } + + kernel( cl::EnqueueArgs( + queue, + cl::NDRange( xSize ) + ), + *scanLineBuffer[rowIndex[0]], + *scanLineBuffer[rowIndex[1]], + *scanLineBuffer[rowIndex[2]], + resultLineBuffer, + rasterParamsBuffer + ); + + queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, inputSize, resultLine.get() ); + + if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine.get(), xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + QgsDebugMsg( "Raster IO Error" ); + } + std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() ); + } + + if ( feedback && feedback->isCanceled() ) + { + //delete the dataset without closing (because it is faster) + gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile ); + return 7; + } + return 0; +} +#endif + + +// TODO: return an anum instead of an int +int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback ) +{ + + GDALAllRegister(); + + //open input file + int xSize, ySize; + gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) ); + if ( !inputDataset ) + { + return 1; //opening of input file failed + } + + //output driver + GDALDriverH outputDriver = openOutputDriver(); + if ( !outputDriver ) + { + return 2; + } + + gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) ); + if ( !outputDataset ) + { + return 3; //create operation on output file failed + } + + //open first raster band for reading (operation is only for single band raster) + GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 ); + if ( !rasterBand ) + { + return 4; + } + mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr ); + + GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 ); + if ( !outputRasterBand ) + { + return 5; + } + //try to set -9999 as nodata value + GDALSetRasterNoDataValue( outputRasterBand, -9999 ); + mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr ); + + if ( ySize < 3 ) //we require at least three rows (should be true for most datasets) + { + return 6; + } + + //keep only three scanlines in memory at a time, make room for initial and final nodata + std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) ); + float *scanLine1 = ( float * ) CPLMalloc( bufferSize ); + float *scanLine2 = ( float * ) CPLMalloc( bufferSize ); + float *scanLine3 = ( float * ) CPLMalloc( bufferSize ); + + 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 yIndex = 0; yIndex < ySize; ++yIndex ) + { + if ( feedback && feedback->isCanceled() ) + { + break; + } + + if ( feedback ) + { + feedback->setProgress( 100.0 * static_cast< double >( yIndex ) / ySize ); + } + + 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 ) + { + scanLine1[a] = mInputNodataValue; + } + // Read scanline2 + if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, &scanLine2[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + QgsDebugMsg( "Raster IO Error" ); + } + } + else + { + //normally fetch only scanLine3 and release scanline 1 if we move forward one row + CPLFree( scanLine1 ); + scanLine1 = scanLine2; + scanLine2 = scanLine3; + scanLine3 = ( float * ) CPLMalloc( bufferSize ); + } + + // Read scanline 3 + if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values + { + for ( int a = 0; a < xSize + 2; ++a ) + { + scanLine3[a] = mInputNodataValue; + } + } + else + { + if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + QgsDebugMsg( "Raster IO Error" ); + } + } + // Set first and last extra columns to nodata + scanLine1[0] = scanLine1[xSize + 1] = mInputNodataValue; + scanLine2[0] = scanLine2[xSize + 1] = mInputNodataValue; + scanLine3[0] = scanLine3[xSize + 1] = mInputNodataValue; + + + + // j is the x axis index, skip 0 and last cell that have been filled with nodata + for ( int xIndex = 0; xIndex < xSize ; ++xIndex ) + { + // 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, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + QgsDebugMsg( "Raster IO Error" ); + } + } + + CPLFree( resultLine ); + CPLFree( scanLine1 ); + CPLFree( scanLine2 ); + CPLFree( scanLine3 ); + + if ( feedback && feedback->isCanceled() ) + { + //delete the dataset without closing (because it is faster) + gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile ); + return 7; + } + return 0; +} diff --git a/src/analysis/raster/qgsninecellfilter.h b/src/analysis/raster/qgsninecellfilter.h index 2442176deb5..63b664708b4 100644 --- a/src/analysis/raster/qgsninecellfilter.h +++ b/src/analysis/raster/qgsninecellfilter.h @@ -22,6 +22,7 @@ #include "gdal.h" #include "qgis_analysis.h" #include "qgsogrutils.h" + class QgsFeedback; /** @@ -39,8 +40,9 @@ class ANALYSIS_EXPORT QgsNineCellFilter /** * Starts the calculation, reads from mInputFile and stores the result in mOutputFile - \param feedback feedback object that receives update and that is checked for cancelation. - \returns 0 in case of success*/ + * \param feedback feedback object that receives update and that is checked for cancelation. + * \returns 0 in case of success + */ int processRaster( QgsFeedback *feedback = nullptr ); double cellSizeX() const { return mCellSizeX; } @@ -57,8 +59,23 @@ class ANALYSIS_EXPORT QgsNineCellFilter void setOutputNodataValue( double value ) { mOutputNodataValue = value; } /** - * 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 + */ virtual float processNineCellWindow( float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) = 0; @@ -80,6 +97,40 @@ class ANALYSIS_EXPORT QgsNineCellFilter \returns the output dataset or nullptr in case of error*/ gdal::dataset_unique_ptr openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver ); + /** + * \brief processRasterCPU executes the computation on the CPU + * \param feedback instance of QgsFeedback, to allow for progress monitoring and cancelation + * \return an opaque integer for error codes: 0 in case of success + */ + int processRasterCPU( QgsFeedback *feedback = nullptr ); + +#ifdef HAVE_OPENCL + + /** + * \brief processRasterGPU executes the computation on the GPU + * \param source path to the OpenCL source file + * \param feedback instance of QgsFeedback, to allow for progress monitoring and cancelation + * \return an opaque integer for error codes: 0 in case of success + */ + int processRasterGPU( const QString &source, QgsFeedback *feedback = nullptr ); + + /** + * \brief addExtraRasterParams allow derived classes to add parameters needed + * by OpenCL program + * \param params vector of parameters passed to OpenCL algorithm + */ + virtual void addExtraRasterParams( std::vector ¶ms ) + { + Q_UNUSED( params ); + } + + virtual const QString openClProgramBaseName() const + { + return QString(); + } + +#endif + protected: QString mInputFile; diff --git a/src/analysis/raster/qgsruggednessfilter.cpp b/src/analysis/raster/qgsruggednessfilter.cpp index 6ef6f7f592a..4440c22f905 100644 --- a/src/analysis/raster/qgsruggednessfilter.cpp +++ b/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 ) diff --git a/src/analysis/raster/qgsruggednessfilter.h b/src/analysis/raster/qgsruggednessfilter.h index 4456a797bce..3dcf5d4c9a1 100644 --- a/src/analysis/raster/qgsruggednessfilter.h +++ b/src/analysis/raster/qgsruggednessfilter.h @@ -39,8 +39,16 @@ class ANALYSIS_EXPORT QgsRuggednessFilter: public QgsNineCellFilter float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) override; +#ifdef HAVE_OPENCL private: QgsRuggednessFilter(); + + virtual const QString openClProgramBaseName() const override + { + return QStringLiteral( "ruggedness" ); + } +#endif + }; #endif // QGSRUGGEDNESSFILTER_H diff --git a/src/analysis/raster/qgsslopefilter.cpp b/src/analysis/raster/qgsslopefilter.cpp index a58fc981700..0b066348e00 100644 --- a/src/analysis/raster/qgsslopefilter.cpp +++ b/src/analysis/raster/qgsslopefilter.cpp @@ -24,8 +24,10 @@ QgsSlopeFilter::QgsSlopeFilter( const QString &inputFile, const QString &outputF } -float QgsSlopeFilter::processNineCellWindow( float *x11, float *x21, float *x31, - float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) +float QgsSlopeFilter::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 ); diff --git a/src/analysis/raster/qgsslopefilter.h b/src/analysis/raster/qgsslopefilter.h index 7c36654971b..d6e843b9d27 100644 --- a/src/analysis/raster/qgsslopefilter.h +++ b/src/analysis/raster/qgsslopefilter.h @@ -35,6 +35,17 @@ class ANALYSIS_EXPORT QgsSlopeFilter: public QgsDerivativeFilter float processNineCellWindow( float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) override; + + +#ifdef HAVE_OPENCL + private: + + virtual const QString openClProgramBaseName() const override + { + return QStringLiteral( "slope" ); + } +#endif + }; #endif // QGSSLOPEFILTER_H diff --git a/src/app/main.cpp b/src/app/main.cpp index 6b9689201ef..49539c41925 100644 --- a/src/app/main.cpp +++ b/src/app/main.cpp @@ -104,6 +104,10 @@ typedef SInt32 SRefCon; #include "qgsuserprofilemanager.h" #include "qgsuserprofile.h" +#ifdef HAVE_OPENCL +#include "qgsopenclutils.h" +#endif + /** * Print usage text */ @@ -141,6 +145,9 @@ void usage( const QString &appName ) << QStringLiteral( "\t[--profile name]\tload a named profile from the users profiles folder.\n" ) << QStringLiteral( "\t[--profiles-path path]\tpath to store user profile folders. Will create profiles inside a {path}\\profiles folder \n" ) << QStringLiteral( "\t[--version-migration]\tforce the settings migration from older version if found\n" ) +#ifdef HAVE_OPENCL + << QStringLiteral( "\t[--openclprogramfolder]\t\tpath to the folder containing the sources for OpenCL programs.\n" ) +#endif << QStringLiteral( "\t[--help]\t\tthis text\n" ) << QStringLiteral( "\t[--]\t\ttreat all following arguments as FILEs\n\n" ) << QStringLiteral( " FILE:\n" ) @@ -543,6 +550,10 @@ int main( int argc, char *argv[] ) QString customizationfile; QString globalsettingsfile; +#ifdef HAVE_OPENCL + QString openClProgramFolder; +#endif + // TODO Fix android #if defined(ANDROID) QgsDebugMsg( QString( "Android: All params stripped" ) );// Param %1" ).arg( argv[0] ) ); @@ -726,6 +737,12 @@ int main( int argc, char *argv[] ) { dxfMapTheme = args[++i]; } +#ifdef HAVE_OPENCL + else if ( arg == QLatin1String( "--openclprogramfolder" ) ) + { + openClProgramFolder = args[++i]; + } +#endif else if ( arg == QLatin1String( "--" ) ) { for ( i++; i < args.size(); ++i ) @@ -1444,6 +1461,24 @@ int main( int argc, char *argv[] ) // make sure we don't have a dirty blank project after launch QgsProject::instance()->setDirty( false ); +#ifdef HAVE_OPENCL + + // Overrides the OpenCL path to the folder containing the programs + // - use the path specified with --openclprogramfolder, + // - use the environment QGIS_OPENCL_PROGRAM_FOLDER if not found + // - use a default location as a fallback (this is set in QgsApplication initialization) + if ( openClProgramFolder.isEmpty() ) + { + openClProgramFolder = getenv( "QGIS_OPENCL_PROGRAM_FOLDER" ); + } + + if ( ! openClProgramFolder.isEmpty() ) + { + QgsOpenClUtils::setSourcePath( openClProgramFolder ); + } + +#endif + ///////////////////////////////////////////////////////////////////// // Continue on to interactive gui... ///////////////////////////////////////////////////////////////////// diff --git a/src/app/qgisapp.cpp b/src/app/qgisapp.cpp index 0042861ebbc..a82e1908a94 100644 --- a/src/app/qgisapp.cpp +++ b/src/app/qgisapp.cpp @@ -97,6 +97,10 @@ #include "qgsgui.h" #include "qgsnative.h" +#ifdef HAVE_OPENCL +#include "qgsopenclutils.h" +#endif + #include #include #include @@ -1103,6 +1107,13 @@ QgisApp::QgisApp( QSplashScreen *splash, bool restorePlugins, bool skipVersionCh // Setup QgsNetworkAccessManager (this needs to happen after authentication, for proxy settings) namSetup(); + +#ifdef HAVE_OPENCL + // Setup the default OpenCL programs source path, this my be overridden later by main.cpp startup + QgsOpenClUtils::setSourcePath( QDir( QgsApplication::pkgDataPath() ).absoluteFilePath( QStringLiteral( "resources/opencl_programs" ) ) ); +#endif + + QgsApplication::dataItemProviderRegistry()->addProvider( new QgsQlrDataItemProvider() ); registerCustomDropHandler( new QgsQlrDropHandler() ); QgsApplication::dataItemProviderRegistry()->addProvider( new QgsQptDataItemProvider() ); diff --git a/src/app/qgsoptions.cpp b/src/app/qgsoptions.cpp index 5a50211e1f4..f4fe9ebe5ab 100644 --- a/src/app/qgsoptions.cpp +++ b/src/app/qgsoptions.cpp @@ -51,6 +51,10 @@ #include "qgslocatorwidget.h" #include "qgslocatoroptionswidget.h" +#ifdef HAVE_OPENCL +#include "qgsopenclutils.h" +#endif + #include #include #include @@ -1076,6 +1080,44 @@ QgsOptions::QgsOptions( QWidget *parent, Qt::WindowFlags fl, const QListaddWidget( page ); } +#ifdef HAVE_OPENCL + + // Setup OpenCL (GPU) widget + mGPUEnableCheckBox->setChecked( QgsOpenClUtils::enabled( ) ); + if ( QgsOpenClUtils::available( ) ) + { + mGPUEnableCheckBox->setEnabled( true ); + + for ( const auto &dev : QgsOpenClUtils::devices( ) ) + { + mOpenClDevicesCombo->addItem( QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Name, dev ), QgsOpenClUtils::deviceId( dev ) ); + } + // Info updater + std::function infoUpdater = [ = ]( int ) + { + mGPUInfoTextBrowser->setText( QgsOpenClUtils::deviceDescription( mOpenClDevicesCombo->currentData().toString() ) ); + }; + connect( mOpenClDevicesCombo, qgis::overload< int >::of( &QComboBox::currentIndexChanged ), infoUpdater ); + mOpenClDevicesCombo->setCurrentIndex( mOpenClDevicesCombo->findData( QgsOpenClUtils::deviceId( QgsOpenClUtils::activeDevice() ) ) ); + infoUpdater( -1 ); + } + else + { + mGPUEnableCheckBox->setEnabled( false ); + mGPUInfoTextBrowser->setText( tr( "An OpenCL compatible device was not found on your system.
" + "You may need to install additional libraries in order to enable OpenCL.
" + "Please check your logs for further details." ) ); + } + + +#else + + mOptionsListWidget->removeItemWidget( mOptionsListWidget->findItems( tr( "Acceleration" ), Qt::MatchExactly ).first() ); + mOptionsStackedWidget->removeWidget( mOptionsPageAcceleration ); + + +#endif + connect( pbnEditCreateOptions, &QAbstractButton::pressed, this, &QgsOptions::editCreateOptions ); connect( pbnEditPyramidsOptions, &QAbstractButton::pressed, this, &QgsOptions::editPyramidsOptions ); @@ -1617,6 +1659,13 @@ void QgsOptions::saveOptions() // Number settings mSettings->setValue( QStringLiteral( "locale/showGroupSeparator" ), cbShowGroupSeparator->isChecked( ) ); +#ifdef HAVE_OPENCL + // OpenCL settings + QgsOpenClUtils::setEnabled( mGPUEnableCheckBox->isChecked() ); + QString preferredDevice( mOpenClDevicesCombo->currentData().toString() ); + QgsOpenClUtils::storePreferredDevice( preferredDevice ); +#endif + // Gdal skip driver list if ( mLoadedGdalDriverList ) saveGdalDriverList(); diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 84146225fc1..399fb80a8a4 100755 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -521,7 +521,6 @@ SET(QGIS_CORE_SRCS qgsuserprofilemanager.cpp ) - FILE(GLOB JSON_HELP_FILES "${CMAKE_SOURCE_DIR}/resources/function_help/json/*") IF(NOT USING_NINJA) STRING(REPLACE "$" "$$" JSON_HELP_FILES "${JSON_HELP_FILES}") @@ -769,8 +768,16 @@ IF (QT_MOBILITY_LOCATION_FOUND OR Qt5Positioning_FOUND) ) ENDIF (QT_MOBILITY_LOCATION_FOUND OR Qt5Positioning_FOUND) + +IF (HAVE_OPENCL) + SET(QGIS_CORE_MOC_HDRS ${QGIS_CORE_MOC_HDRS} + qgsopenclutils.h + ) +ENDIF (HAVE_OPENCL) + QT5_WRAP_CPP(QGIS_CORE_MOC_SRCS ${QGIS_CORE_MOC_HDRS}) + IF(MSVC) SET_SOURCE_FILES_PROPERTIES(${QGIS_CORE_MOC_SRCS} PROPERTIES COMPILE_FLAGS "/wd4512 /wd4996" ) ELSE(MSVC) @@ -1210,6 +1217,17 @@ INCLUDE_DIRECTORIES(SYSTEM ) +IF (HAVE_OPENCL) + SET(QGIS_CORE_SRCS ${QGIS_CORE_SRCS} + qgsopenclutils.cpp + ) + SET(QGIS_CORE_HDRS ${QGIS_CORE_HDRS} + qgsopenclutils.h + ) + INCLUDE_DIRECTORIES(${OpenCL_INCLUDE_DIR}) +ENDIF (HAVE_OPENCL) + + IF (APPLE) # Libtasn1 is for DER-encoded PKI ASN.1 parsing/extracting workarounds INCLUDE_DIRECTORIES(SYSTEM @@ -1245,6 +1263,12 @@ IF(ENABLE_MODELTEST) TARGET_LINK_LIBRARIES(qgis_core ${Qt5Test_LIBRARIES}) ENDIF(ENABLE_MODELTEST) + +IF(HAVE_OPENCL) + TARGET_LINK_LIBRARIES(qgis_core ${OpenCL_LIBRARY}) +ENDIF(HAVE_OPENCL) + + IF(NOT APPLE) INSTALL(FILES ${QGIS_CORE_HDRS} ${QGIS_CORE_MOC_HDRS} DESTINATION ${QGIS_INCLUDE_DIR}) ELSE(NOT APPLE) diff --git a/src/core/qgsopenclutils.cpp b/src/core/qgsopenclutils.cpp new file mode 100644 index 00000000000..fae95bcac12 --- /dev/null +++ b/src/core/qgsopenclutils.cpp @@ -0,0 +1,530 @@ +/*************************************************************************** + qgsopenclutils.cpp - QgsOpenClUtils + + --------------------- + begin : 11.4.2018 + copyright : (C) 2018 by elpaso + email : elpaso at itopen dot it + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +#include "qgsopenclutils.h" +#include "qgssettings.h" +#include "qgsmessagelog.h" +#include "qgslogger.h" + + +#include +#include +#include + +QLatin1String QgsOpenClUtils::SETTINGS_GLOBAL_ENABLED_KEY = QLatin1Literal( "OpenClEnabled" ); +QLatin1String QgsOpenClUtils::SETTINGS_DEFAULT_DEVICE_KEY = QLatin1Literal( "OpenClDefaultDevice" ); +QLatin1String QgsOpenClUtils::LOGMESSAGE_TAG = QLatin1Literal( "OpenCL" ); +bool QgsOpenClUtils::sAvailable = false; +cl::Platform QgsOpenClUtils::sDefaultPlatform = cl::Platform(); +cl::Device QgsOpenClUtils::sActiveDevice = cl::Device(); +QString QgsOpenClUtils::sSourcePath = QString(); + + +const std::vector QgsOpenClUtils::devices() +{ + std::vector platforms; + cl::Platform::get( &platforms ); + std::vector existingDevices; + for ( auto &p : platforms ) + { + std::string platver = p.getInfo(); + QgsDebugMsg( QStringLiteral( "Found OpenCL platform %1: %2" ).arg( QString::fromStdString( platver ), QString::fromStdString( p.getInfo() ) ) ); + if ( platver.find( "OpenCL 1." ) != std::string::npos ) + { + std::vector _devices; + // Check for a device + try + { + p.getDevices( CL_DEVICE_TYPE_ALL, &_devices ); + } + catch ( cl::Error &e ) + { + QgsDebugMsgLevel( QStringLiteral( "Error %1 on platform %3 searching for OpenCL device: %2" ) + .arg( errorText( e.err() ), + QString::fromStdString( e.what() ), + QString::fromStdString( p.getInfo() ) ), 2 ); + } + if ( _devices.size() > 0 ) + { + for ( unsigned long i = 0; i < _devices.size(); i++ ) + { + existingDevices.push_back( _devices[i] ); + } + } + } + } + return existingDevices; +} + +void QgsOpenClUtils::init() +{ + static std::once_flag initialized; + std::call_once( initialized, [ = ]( ) + { + try + { + activate( preferredDevice() ); + } + catch ( cl::Error &e ) + { + QgsMessageLog::logMessage( QObject::tr( "Error %1 initializing OpenCL device: %2" ) + .arg( errorText( e.err() ), QString::fromStdString( e.what() ) ), + LOGMESSAGE_TAG, Qgis::Critical ); + } + + } ); +} + +QString QgsOpenClUtils::sourcePath() +{ + return sSourcePath; +} + +void QgsOpenClUtils::setSourcePath( const QString &value ) +{ + sSourcePath = value; +} + +QString QgsOpenClUtils::activeDeviceInfo( const QgsOpenClUtils::Info infoType ) +{ + return deviceInfo( infoType, activeDevice( ) ); +} + +QString QgsOpenClUtils::deviceInfo( const Info infoType, cl::Device device ) +{ + try + { + switch ( infoType ) + { + case Info::Vendor: + return QString::fromStdString( device.getInfo() ); + case Info::Profile: + return QString::fromStdString( device.getInfo() ); + case Info::Version: + return QString::fromStdString( device.getInfo() ); + case Info::ImageSupport: + return device.getInfo() ? QStringLiteral( "True" ) : QStringLiteral( "False" ); + case Info::Image2dMaxHeight: + return QString::number( device.getInfo() ); + case Info::MaxMemAllocSize: + return QString::number( device.getInfo() ); + case Info::Image2dMaxWidth: + return QString::number( device.getInfo() ); + case Info::Type: + { + unsigned long type( device.getInfo() ); + int mappedType; + switch ( type ) + { + case CL_DEVICE_TYPE_CPU: + mappedType = QgsOpenClUtils::HardwareType::CPU; + break; + case CL_DEVICE_TYPE_GPU: + mappedType = QgsOpenClUtils::HardwareType::GPU; + break; + default: + mappedType = QgsOpenClUtils::HardwareType::Other; + } + QMetaEnum metaEnum = QMetaEnum::fromType(); + return metaEnum.valueToKey( mappedType ); + } + case Info::Name: + return QString::fromStdString( device.getInfo() ); + } + } + catch ( cl::Error &e ) + { + // This can be a legitimate error when initializing, let's log it quietly + QgsDebugMsgLevel( QStringLiteral( "Error %1 getting info for OpenCL device: %2" ) + .arg( errorText( e.err() ), QString::fromStdString( e.what() ) ), + 4 ); + return QString(); + } +} + + +bool QgsOpenClUtils::enabled() +{ + return QgsSettings().value( SETTINGS_GLOBAL_ENABLED_KEY, false, QgsSettings::Section::Core ).toBool(); +} + +cl::Device QgsOpenClUtils::activeDevice() +{ + return sActiveDevice; +} + +void QgsOpenClUtils::storePreferredDevice( const QString deviceId ) +{ + QgsSettings().setValue( SETTINGS_DEFAULT_DEVICE_KEY, deviceId, QgsSettings::Section::Core ); +} + +QString QgsOpenClUtils::preferredDevice() +{ + return QgsSettings().value( SETTINGS_DEFAULT_DEVICE_KEY, QString( ), QgsSettings::Section::Core ).toString(); +} + +QString QgsOpenClUtils::deviceId( const cl::Device device ) +{ + return QStringLiteral( "%1|%2|%3|%4" ) + .arg( deviceInfo( QgsOpenClUtils::Info::Name, device ) ) + .arg( deviceInfo( QgsOpenClUtils::Info::Vendor, device ) ) + .arg( deviceInfo( QgsOpenClUtils::Info::Version, device ) ) + .arg( deviceInfo( QgsOpenClUtils::Info::Type, device ) ); +} + +bool QgsOpenClUtils::activate( const QString &preferredDeviceId ) +{ + if ( deviceId( activeDevice() ) == preferredDeviceId ) + { + return false; + } + try + { + std::vector platforms; + cl::Platform::get( &platforms ); + cl::Platform plat; + cl::Device dev; + bool deviceFound = false; + for ( auto &p : platforms ) + { + if ( deviceFound ) + break; + std::string platver = p.getInfo(); + QgsDebugMsg( QStringLiteral( "Found OpenCL platform %1: %2" ).arg( QString::fromStdString( platver ), QString::fromStdString( p.getInfo() ) ) ); + if ( platver.find( "OpenCL 1." ) != std::string::npos ) + { + std::vector devices; + // Search for a device + try + { + p.getDevices( CL_DEVICE_TYPE_ALL, &devices ); + // First search for the preferred device + if ( ! preferredDeviceId.isEmpty() ) + { + for ( const auto &_dev : devices ) + { + if ( preferredDeviceId == deviceId( _dev ) ) + { + // Got one! + plat = p; + dev = _dev; + deviceFound = true; + break; + } + } + } + // Not found or preferred device id not set: get the first GPU + if ( ! deviceFound ) + { + for ( const auto &_dev : devices ) + { + if ( preferredDeviceId.isEmpty() && _dev.getInfo() == CL_DEVICE_TYPE_GPU ) + { + // Got one! + plat = p; + dev = _dev; + deviceFound = true; + break; + } + } + } + // Still nothing? Get the first device + if ( ! deviceFound ) + { + for ( const auto &_dev : devices ) + { + if ( preferredDeviceId.isEmpty() && _dev.getInfo() == CL_DEVICE_TYPE_GPU ) + { + // Got one! + plat = p; + dev = _dev; + deviceFound = true; + break; + } + } + } + if ( ! deviceFound ) + { + QgsMessageLog::logMessage( QObject::tr( "No OpenCL 1.x device could be found." ), LOGMESSAGE_TAG, Qgis::Warning ); + } + } + catch ( cl::Error &e ) + { + QgsDebugMsg( QStringLiteral( "Error %1 on platform %3 searching for OpenCL device: %2" ) + .arg( errorText( e.err() ), + QString::fromStdString( e.what() ), + QString::fromStdString( p.getInfo() ) ) ); + } + + } + } + if ( ! plat() ) + { + QgsMessageLog::logMessage( QObject::tr( "No OpenCL 1.x platform found." ), LOGMESSAGE_TAG, Qgis::Warning ); + sAvailable = false; + } + else + { + cl::Platform newP = cl::Platform::setDefault( plat ); + if ( newP != plat ) + { + QgsMessageLog::logMessage( QObject::tr( "Error setting default platform." ), + LOGMESSAGE_TAG, Qgis::Warning ); + sAvailable = false; + } + else + { + cl::Device::setDefault( dev ); + QgsMessageLog::logMessage( QObject::tr( "Active OpenCL device: %1" ) + .arg( QString::fromStdString( dev.getInfo() ) ), + LOGMESSAGE_TAG, Qgis::Success ); + sAvailable = true; + sActiveDevice = dev; + sDefaultPlatform = plat; + } + } + } + + catch ( cl::Error &e ) + { + QgsMessageLog::logMessage( QObject::tr( "Error %1 searching for OpenCL device: %2" ) + .arg( errorText( e.err() ), QString::fromStdString( e.what() ) ), + LOGMESSAGE_TAG, Qgis::Critical ); + sAvailable = false; + } + return sAvailable; +} + +QString QgsOpenClUtils::deviceDescription( const cl::Device device ) +{ + return QStringLiteral( + "Type: %9
" + "Name: %1
" + "Vendor: %2
" + "Profile: %3
" + "Version: %4
" + "Image support: %5
" + "Max image2d width: %6
" + "Max image2d height: %7
" + "Max mem alloc size: %8
" + ).arg( QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Name, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Vendor, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Profile, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Version, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::ImageSupport, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Image2dMaxWidth, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Image2dMaxHeight, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::MaxMemAllocSize, device ), + QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Type, device ) ); +} + +QString QgsOpenClUtils::deviceDescription( const QString deviceId ) +{ + for ( const auto &dev : devices( ) ) + { + if ( QgsOpenClUtils::deviceId( dev ) == deviceId ) + return deviceDescription( dev ); + } + return QString(); +} + +bool QgsOpenClUtils::available() +{ + init(); + return sAvailable; +} + +void QgsOpenClUtils::setEnabled( bool enabled ) +{ + QgsSettings().setValue( SETTINGS_GLOBAL_ENABLED_KEY, enabled, QgsSettings::Section::Core ); +} + + + +QString QgsOpenClUtils::sourceFromPath( const QString &path ) +{ + // TODO: check for compatibility with current platform ( cl_khr_fp64 ) + // Try to load the program sources + QString source_str; + QFile file( path ); + if ( file.open( QFile::ReadOnly | QFile::Text ) ) + { + QTextStream in( &file ); + source_str = in.readAll(); + file.close(); + } + else + { + QgsMessageLog::logMessage( QObject::tr( "Could not load OpenCL program from path %1." ).arg( path ), LOGMESSAGE_TAG, Qgis::Warning ); + } + return source_str; +} + +QString QgsOpenClUtils::sourceFromBaseName( const QString &baseName ) +{ + QString path = QStringLiteral( "%1/%2.cl" ).arg( sourcePath(), baseName ); + return sourceFromPath( path ); +} + +QString QgsOpenClUtils::buildLog( cl::BuildError &error ) +{ + cl::BuildLogType build_logs = error.getBuildLog(); + QString build_log; + if ( build_logs.size() > 0 ) + build_log = QString::fromStdString( build_logs[0].second ); + return build_log; +} + +QString QgsOpenClUtils::errorText( const int errorCode ) +{ + switch ( errorCode ) + { + case 0: return QStringLiteral( "CL_SUCCESS" ); + case -1: return QStringLiteral( "CL_DEVICE_NOT_FOUND" ); + case -2: return QStringLiteral( "CL_DEVICE_NOT_AVAILABLE" ); + case -3: return QStringLiteral( "CL_COMPILER_NOT_AVAILABLE" ); + case -4: return QStringLiteral( "CL_MEM_OBJECT_ALLOCATION_FAILURE" ); + case -5: return QStringLiteral( "CL_OUT_OF_RESOURCES" ); + case -6: return QStringLiteral( "CL_OUT_OF_HOST_MEMORY" ); + case -7: return QStringLiteral( "CL_PROFILING_INFO_NOT_AVAILABLE" ); + case -8: return QStringLiteral( "CL_MEM_COPY_OVERLAP" ); + case -9: return QStringLiteral( "CL_IMAGE_FORMAT_MISMATCH" ); + case -10: return QStringLiteral( "CL_IMAGE_FORMAT_NOT_SUPPORTED" ); + case -12: return QStringLiteral( "CL_MAP_FAILURE" ); + case -13: return QStringLiteral( "CL_MISALIGNED_SUB_BUFFER_OFFSET" ); + case -14: return QStringLiteral( "CL_EXEC_STATUS_ERROR_FOR_EVENTS_IN_WAIT_LIST" ); + case -15: return QStringLiteral( "CL_COMPILE_PROGRAM_FAILURE" ); + case -16: return QStringLiteral( "CL_LINKER_NOT_AVAILABLE" ); + case -17: return QStringLiteral( "CL_LINK_PROGRAM_FAILURE" ); + case -18: return QStringLiteral( "CL_DEVICE_PARTITION_FAILED" ); + case -19: return QStringLiteral( "CL_KERNEL_ARG_INFO_NOT_AVAILABLE" ); + case -30: return QStringLiteral( "CL_INVALID_VALUE" ); + case -31: return QStringLiteral( "CL_INVALID_DEVICE_TYPE" ); + case -32: return QStringLiteral( "CL_INVALID_PLATFORM" ); + case -33: return QStringLiteral( "CL_INVALID_DEVICE" ); + case -34: return QStringLiteral( "CL_INVALID_CONTEXT" ); + case -35: return QStringLiteral( "CL_INVALID_QUEUE_PROPERTIES" ); + case -36: return QStringLiteral( "CL_INVALID_COMMAND_QUEUE" ); + case -37: return QStringLiteral( "CL_INVALID_HOST_PTR" ); + case -38: return QStringLiteral( "CL_INVALID_MEM_OBJECT" ); + case -39: return QStringLiteral( "CL_INVALID_IMAGE_FORMAT_DESCRIPTOR" ); + case -40: return QStringLiteral( "CL_INVALID_IMAGE_SIZE" ); + case -41: return QStringLiteral( "CL_INVALID_SAMPLER" ); + case -42: return QStringLiteral( "CL_INVALID_BINARY" ); + case -43: return QStringLiteral( "CL_INVALID_BUILD_OPTIONS" ); + case -44: return QStringLiteral( "CL_INVALID_PROGRAM" ); + case -45: return QStringLiteral( "CL_INVALID_PROGRAM_EXECUTABLE" ); + case -46: return QStringLiteral( "CL_INVALID_KERNEL_NAME" ); + case -47: return QStringLiteral( "CL_INVALID_KERNEL_DEFINITION" ); + case -48: return QStringLiteral( "CL_INVALID_KERNEL" ); + case -49: return QStringLiteral( "CL_INVALID_ARG_INDEX" ); + case -50: return QStringLiteral( "CL_INVALID_ARG_VALUE" ); + case -51: return QStringLiteral( "CL_INVALID_ARG_SIZE" ); + case -52: return QStringLiteral( "CL_INVALID_KERNEL_ARGS" ); + case -53: return QStringLiteral( "CL_INVALID_WORK_DIMENSION" ); + case -54: return QStringLiteral( "CL_INVALID_WORK_GROUP_SIZE" ); + case -55: return QStringLiteral( "CL_INVALID_WORK_ITEM_SIZE" ); + case -56: return QStringLiteral( "CL_INVALID_GLOBAL_OFFSET" ); + case -57: return QStringLiteral( "CL_INVALID_EVENT_WAIT_LIST" ); + case -58: return QStringLiteral( "CL_INVALID_EVENT" ); + case -59: return QStringLiteral( "CL_INVALID_OPERATION" ); + case -60: return QStringLiteral( "CL_INVALID_GL_OBJECT" ); + case -61: return QStringLiteral( "CL_INVALID_BUFFER_SIZE" ); + case -62: return QStringLiteral( "CL_INVALID_MIP_LEVEL" ); + case -63: return QStringLiteral( "CL_INVALID_GLOBAL_WORK_SIZE" ); + case -64: return QStringLiteral( "CL_INVALID_PROPERTY" ); + case -65: return QStringLiteral( "CL_INVALID_IMAGE_DESCRIPTOR" ); + case -66: return QStringLiteral( "CL_INVALID_COMPILER_OPTIONS" ); + case -67: return QStringLiteral( "CL_INVALID_LINKER_OPTIONS" ); + case -68: return QStringLiteral( "CL_INVALID_DEVICE_PARTITION_COUNT" ); + case -69: return QStringLiteral( "CL_INVALID_PIPE_SIZE" ); + case -70: return QStringLiteral( "CL_INVALID_DEVICE_QUEUE" ); + case -71: return QStringLiteral( "CL_INVALID_SPEC_ID" ); + case -72: return QStringLiteral( "CL_MAX_SIZE_RESTRICTION_EXCEEDED" ); + case -1002: return QStringLiteral( "CL_INVALID_D3D10_DEVICE_KHR" ); + case -1003: return QStringLiteral( "CL_INVALID_D3D10_RESOURCE_KHR" ); + case -1004: return QStringLiteral( "CL_D3D10_RESOURCE_ALREADY_ACQUIRED_KHR" ); + case -1005: return QStringLiteral( "CL_D3D10_RESOURCE_NOT_ACQUIRED_KHR" ); + case -1006: return QStringLiteral( "CL_INVALID_D3D11_DEVICE_KHR" ); + case -1007: return QStringLiteral( "CL_INVALID_D3D11_RESOURCE_KHR" ); + case -1008: return QStringLiteral( "CL_D3D11_RESOURCE_ALREADY_ACQUIRED_KHR" ); + case -1009: return QStringLiteral( "CL_D3D11_RESOURCE_NOT_ACQUIRED_KHR" ); + case -1010: return QStringLiteral( "CL_INVALID_DX9_MEDIA_ADAPTER_KHR" ); + case -1011: return QStringLiteral( "CL_INVALID_DX9_MEDIA_SURFACE_KHR" ); + case -1012: return QStringLiteral( "CL_DX9_MEDIA_SURFACE_ALREADY_ACQUIRED_KHR" ); + case -1013: return QStringLiteral( "CL_DX9_MEDIA_SURFACE_NOT_ACQUIRED_KHR" ); + case -1093: return QStringLiteral( "CL_INVALID_EGL_OBJECT_KHR" ); + case -1092: return QStringLiteral( "CL_EGL_RESOURCE_NOT_ACQUIRED_KHR" ); + case -1001: return QStringLiteral( "CL_PLATFORM_NOT_FOUND_KHR" ); + case -1057: return QStringLiteral( "CL_DEVICE_PARTITION_FAILED_EXT" ); + case -1058: return QStringLiteral( "CL_INVALID_PARTITION_COUNT_EXT" ); + case -1059: return QStringLiteral( "CL_INVALID_PARTITION_NAME_EXT" ); + case -1094: return QStringLiteral( "CL_INVALID_ACCELERATOR_INTEL" ); + case -1095: return QStringLiteral( "CL_INVALID_ACCELERATOR_TYPE_INTEL" ); + case -1096: return QStringLiteral( "CL_INVALID_ACCELERATOR_DESCRIPTOR_INTEL" ); + case -1097: return QStringLiteral( "CL_ACCELERATOR_TYPE_NOT_SUPPORTED_INTEL" ); + case -1000: return QStringLiteral( "CL_INVALID_GL_SHAREGROUP_REFERENCE_KHR" ); + case -1098: return QStringLiteral( "CL_INVALID_VA_API_MEDIA_ADAPTER_INTEL" ); + case -1099: return QStringLiteral( "CL_INVALID_VA_API_MEDIA_SURFACE_INTEL" ); + case -1100: return QStringLiteral( "CL_VA_API_MEDIA_SURFACE_ALREADY_ACQUIRED_INTEL" ); + case -1101: return QStringLiteral( "CL_VA_API_MEDIA_SURFACE_NOT_ACQUIRED_INTEL" ); + default: return QStringLiteral( "CL_UNKNOWN_ERROR" ); + } +} + +cl::Context QgsOpenClUtils::context() +{ + static cl::Context context; + static std::once_flag contextCreated; + std::call_once( contextCreated, [ = ]() + { + if ( available() && sDefaultPlatform() && sActiveDevice() ) + { + context = cl::Context( sActiveDevice ); + } + } ); + return context; +} + +cl::Program QgsOpenClUtils::buildProgram( const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior ) +{ + cl::Program program; + try + { + program = cl::Program( context, source.toStdString( ) ); + // OpenCL 1.1 for compatibility with older hardware + // TODO: make this configurable + program.build( QStringLiteral( "-cl-std=CL1.1 -I%1" ).arg( sourcePath() ).toStdString().c_str() ); + } + catch ( cl::BuildError &e ) + { + QString build_log( buildLog( e ) ); + if ( build_log.isEmpty() ) + build_log = QObject::tr( "Build logs not available!" ); + QString err = QObject::tr( "Error building OpenCL program: %1" ) + .arg( build_log ); + QgsMessageLog::logMessage( err, LOGMESSAGE_TAG, Qgis::Critical ); + if ( exceptionBehavior == Throw ) + throw e; + } + catch ( cl::Error &e ) + { + QString err = QObject::tr( "Error %1 building OpenCL program in %2" ) + .arg( errorText( e.err() ), QString::fromStdString( e.what() ) ); + QgsMessageLog::logMessage( err, LOGMESSAGE_TAG, Qgis::Critical ); + throw e; + } + return program; +} diff --git a/src/core/qgsopenclutils.h b/src/core/qgsopenclutils.h new file mode 100644 index 00000000000..b35926d1a48 --- /dev/null +++ b/src/core/qgsopenclutils.h @@ -0,0 +1,289 @@ +/*************************************************************************** + qgsopenclutils.h - QgsOpenClUtils + + --------------------- + begin : 11.4.2018 + copyright : (C) 2018 by Alessandro Pasotti + email : elpaso at itopen dot it + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +#ifndef QGSOPENCLUTILS_H +#define QGSOPENCLUTILS_H + +#define SIP_NO_FILE + +#define CL_HPP_ENABLE_EXCEPTIONS +#define CL_HPP_MINIMUM_OPENCL_VERSION 110 +#define CL_HPP_TARGET_OPENCL_VERSION 110 +#include + +#include "qgis_core.h" +#include "qgis.h" + +#include "cpl_conv.h" + +/** + * \ingroup core + * \class QgsOpenClUtils + * \brief The QgsOpenClUtils class is responsible for common OpenCL operations such as + * - enable/disable opencl + * - store and retrieve preferences for the default device + * - check opencl device availability and automatically choose the first GPU device + * - creating the default context + * - loading program sources from standard locations + * - build programs and log errors + * + * Usage: + * + * \code{.cpp} + // This will check if OpenCL is enabled in user options and if there is a suitable + // device, if a device is found it is initialized. + if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() ) + { + // Use the default context + cl::Context ctx = QgsOpenClUtils::context(); + cl::CommandQueue queue( ctx ); + // Load the program from a standard location and build it + cl::Program program = QgsOpenClUtils::buildProgram( ctx, QgsOpenClUtils::sourceFromBaseName( QStringLiteral ( "hillshade" ) ) ); + // Continue with the usual OpenCL buffer, kernel and execution + ... + } + * \endcode + * + * \note not available in Python bindings + * \since QGIS 3.4 + */ +class CORE_EXPORT QgsOpenClUtils +{ + Q_GADGET + + public: + + /** + * The ExceptionBehavior enum define how exceptions generated by OpenCL should be treated + */ + enum ExceptionBehavior + { + Catch, //! Write errors in the message log and silently fail + Throw //! Write errors in the message log and re-throw exceptions + }; + + /** + * The Type enum represent OpenCL device type + */ + enum HardwareType + { + CPU, + GPU, + Other + }; + + Q_ENUM( HardwareType ) + + /** + * The Info enum maps to OpenCL info constants + * + * \see deviceInfo() + */ + enum Info + { + Name = CL_DEVICE_NAME, + Vendor = CL_DEVICE_VENDOR, + Version = CL_DEVICE_VERSION, + Profile = CL_DEVICE_PROFILE, + ImageSupport = CL_DEVICE_IMAGE_SUPPORT, + Image2dMaxWidth = CL_DEVICE_IMAGE2D_MAX_WIDTH, + Image2dMaxHeight = CL_DEVICE_IMAGE2D_MAX_HEIGHT, + MaxMemAllocSize = CL_DEVICE_MAX_MEM_ALLOC_SIZE, + Type = CL_DEVICE_TYPE // CPU/GPU etc. + }; + + /** + * Checks whether a suitable OpenCL platform and device is available on this system + * and initialize the QGIS OpenCL system by activating the preferred device + * if specified in the user the settings, if no preferred device was set or + * the preferred device could not be found the first GPU device is activated, + * the first CPU device acts as a fallback if none of the previous could be found. + * + * This function must always be called before using QGIS OpenCL utils + */ + static bool available(); + + //! Returns true if OpenCL is enabled in the user settings + static bool enabled(); + + //! Returns a list of OpenCL devices found on this sysytem + static const std::vector devices(); + + /** + * Returns the active device. + * + * The active device is set as the default device for all OpenCL operations, + * once it is set it cannot be changed until QGIS is restarted (this is + * due to the way the underlying OpenCL library is built). + */ + static cl::Device activeDevice( ); + + //! Store in the settings the preferred \a deviceId device identifier + static void storePreferredDevice( const QString deviceId ); + + //! Read from the settings the preferred device identifier + static QString preferredDevice( ); + + //! Create a string identifier from a \a device + static QString deviceId( const cl::Device device ); + + /** + * Returns a formatted description for the \a device + */ + static QString deviceDescription( const cl::Device device ); + + /** + * Returns a formatted description for the device identified by \a deviceId + */ + static QString deviceDescription( const QString deviceId ); + + //! Set the OpenCL user setting to \a enabled + static void setEnabled( bool enabled ); + + //! Extract and return the build log error from \a error + static QString buildLog( cl::BuildError &error ); + + //! Read an OpenCL source file from \a path + static QString sourceFromPath( const QString &path ); + + //! Returns the full path to a an OpenCL source file from the \a baseName ('.cl' extension is automatically appended) + static QString sourceFromBaseName( const QString &baseName ); + + //! OpenCL string for message logs + static QLatin1String LOGMESSAGE_TAG; + + //! Returns a string representation from an OpenCL \a errorCode + static QString errorText( const int errorCode ); + + /** + * Build the program from \a source in the given \a context and depending on \a exceptionBehavior + * can throw or catch the exceptions + * \return the built program + */ + static cl::Program buildProgram( const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior = Catch ); + + /** + * Context factory + * + * \return a new context for the default device or an invalid context if + * no device were identified or OpenCL support is not available + * and enabled + */ + static cl::Context context(); + + //! Returns the base path to OpenCL program directory + static QString sourcePath(); + + //! Set the base path to OpenCL program directory + static void setSourcePath( const QString &value ); + + //! Returns \a infoType information about the active (default) device + static QString activeDeviceInfo( const Info infoType = Info::Name ); + + //! Returns \a infoType information about the \a device + static QString deviceInfo( const Info infoType, cl::Device device ); + + /** + * Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because + * OpenCL C++ API may throw exceptions + */ + template + struct CPLAllocator + { + + public: + + explicit CPLAllocator( unsigned long size ): mMem( static_cast( CPLMalloc( sizeof( T ) * size ) ) ) { } + + ~CPLAllocator() + { + CPLFree( static_cast( mMem ) ); + } + + void reset( T *newData ) + { + if ( mMem ) + CPLFree( static_cast( mMem ) ); + mMem = newData; + } + + void reset( unsigned long size ) + { + reset( static_cast( CPLMalloc( sizeof( T ) *size ) ) ); + } + + T &operator* () + { + return &mMem[0]; + } + + T *release() + { + T *tmpMem = mMem; + mMem = nullptr; + return tmpMem; + } + + T &operator[]( const int index ) + { + return mMem[index]; + } + + T *get() + { + return mMem; + } + + private: + + T *mMem = nullptr; + }; + + + private: + + QgsOpenClUtils(); + + /** + * Activate a device identified by its \a preferredDeviceId by making it the default device + * if the device does not exists or deviceId is empty, the first GPU device will be + * activated, if a GPU device is not found, the first CPU device will be chosen instead. + * + * Called once by init() when OpenCL is used for the first time in a QGIS working session. + * + * \return true if the device could be found and activated. Return false if the device was already + * the active one or if a device could not be activated. + * + * \see init() + * \see available() + */ + static bool activate( const QString &preferredDeviceId = QString() ); + + /** + * Initialize the OpenCL system by setting and activating the default device. + */ + static void init(); + + static bool sAvailable; + static cl::Device sActiveDevice; + static cl::Platform sDefaultPlatform; + static QLatin1String SETTINGS_GLOBAL_ENABLED_KEY; + static QLatin1String SETTINGS_DEFAULT_DEVICE_KEY; + static QString sSourcePath; +}; + + + +#endif // QGSOPENCLUTILS_H diff --git a/src/core/raster/qgscolorrampshader.cpp b/src/core/raster/qgscolorrampshader.cpp index 0148dbaf317..841c71c2728 100644 --- a/src/core/raster/qgscolorrampshader.cpp +++ b/src/core/raster/qgscolorrampshader.cpp @@ -359,6 +359,9 @@ bool QgsColorRampShader::shade( double value, int *returnRedValue, int *returnGr mLUTInitialized = true; } + if ( mLUT.empty() ) + return false; + // overflow indicates that value > maximum value + DOUBLE_DIFF_THRESHOLD // that way idx can point to the last valid item bool overflow = false; @@ -377,10 +380,6 @@ bool QgsColorRampShader::shade( double value, int *returnRedValue, int *returnGr overflow = true; } } - else if ( lutIndex < 0 ) - { - return false; - } else { // get initial value from LUT diff --git a/src/core/raster/qgshillshaderenderer.cpp b/src/core/raster/qgshillshaderenderer.cpp index 38f3d040010..6c0a0bfa937 100644 --- a/src/core/raster/qgshillshaderenderer.cpp +++ b/src/core/raster/qgshillshaderenderer.cpp @@ -22,12 +22,15 @@ #include "qgsrasterinterface.h" #include "qgsrasterblock.h" #include "qgsrectangle.h" -#include "qgssettings.h" +#include "qgsmessagelog.h" #include +#ifdef HAVE_OPENCL #ifdef QGISDEBUG -#include "qgsmessagelog.h" #include +#include "qgssettings.h" +#endif +#include "qgsopenclutils.h" #endif QgsHillshadeRenderer::QgsHillshadeRenderer( QgsRasterInterface *input, int band, double lightAzimuth, double lightAngle ): @@ -128,192 +131,377 @@ QgsRasterBlock *QgsHillshadeRenderer::block( int bandNo, const QgsRectangle &ext return outputBlock.release(); } - float cellXSize = extent.width() / static_cast( width ); - float cellYSize = extent.height() / static_cast( height ); - float zenithRad = std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0; - float azimuthRad = -1 * mLightAzimuth * M_PI / 180.0; + // Starting the computation + + // Common pre-calculated values + float cellXSize = static_cast( extent.width() ) / width; + float cellYSize = static_cast( extent.height() ) / height; + float zenithRad = static_cast( std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0 ); + float azimuthRad = static_cast( -1 * mLightAzimuth * M_PI / 180.0 ); float cosZenithRad = std::cos( zenithRad ); float sinZenithRad = std::sin( zenithRad ); // For fast formula from GDAL DEM - float cos_alt_mul_z = cosZenithRad * mZFactor; + float cos_alt_mul_z = cosZenithRad * static_cast( mZFactor ); float cos_az_mul_cos_alt_mul_z = std::cos( azimuthRad ) * cos_alt_mul_z; float sin_az_mul_cos_alt_mul_z = std::sin( azimuthRad ) * cos_alt_mul_z; - float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0 * cos_az_mul_cos_alt_mul_z; - float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0 * sin_az_mul_cos_alt_mul_z; - float square_z = mZFactor * mZFactor; - float sinZenithRad_mul_254 = 254.0 * sinZenithRad; + float cos_az_mul_cos_alt_mul_z_mul_254 = 254.0f * cos_az_mul_cos_alt_mul_z; + float sin_az_mul_cos_alt_mul_z_mul_254 = 254.0f * sin_az_mul_cos_alt_mul_z; + float square_z = static_cast( mZFactor * mZFactor ); + float sin_altRadians_mul_254 = 254.0f * sinZenithRad; // For multi directional - float sinZenithRad_mul_127 = 127.0 * sinZenithRad; + float sin_altRadians_mul_127 = 127.0f * sinZenithRad; // 127.0 * std::cos(225.0 * M_PI / 180.0) = -32.87001872802012 float cos225_az_mul_cos_alt_mul_z_mul_127 = -32.87001872802012f * cos_alt_mul_z; - float cos_alt_mul_z_mul_127 = 127.0 * cos_alt_mul_z; + float cos_alt_mul_z_mul_127 = 127.0f * cos_alt_mul_z; - QRgb myDefaultColor = NODATA_COLOR; + QRgb defaultNodataColor = NODATA_COLOR; + + +#ifdef HAVE_OPENCL + + // Use OpenCL? For now OpenCL is enabled in the default configuration only + bool useOpenCL( QgsOpenClUtils::enabled() + && QgsOpenClUtils::available() + && ( ! mRasterTransparency || mRasterTransparency->isEmpty() ) + && mAlphaBand <= 0 ); + // Check for sources + QString source; + if ( useOpenCL ) + { + source = QgsOpenClUtils::sourceFromBaseName( QStringLiteral( "hillshade_renderer" ) ); + if ( source.isEmpty() ) + { + useOpenCL = false; + QgsMessageLog::logMessage( QObject::tr( "Error loading OpenCL program source from path" ).arg( QgsOpenClUtils::sourcePath() ), QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical ); + } + } #ifdef QGISDEBUG std::chrono::time_point startTime( std::chrono::system_clock::now() ); #endif - for ( qgssize i = 0; i < ( qgssize )height; i++ ) + if ( useOpenCL ) { - for ( qgssize j = 0; j < ( qgssize )width; j++ ) + try { + std::size_t inputDataTypeSize = inputBlock->dataTypeSize(); + std::size_t outputDataTypeSize = outputBlock->dataTypeSize(); + // Buffer scanline, 1px height, 2px wider + // Data type for input is Float32 (4 bytes) + std::size_t scanLineWidth( inputBlock->width() + 2 ); + std::size_t inputSize( inputDataTypeSize * inputBlock->width() ); - if ( inputBlock->isNoData( i, j ) ) + // CL buffers are also 2px wider for nodata initial and final columns + std::size_t bufferWidth( width + 2 ); + std::size_t bufferSize( inputDataTypeSize * bufferWidth ); + + // Buffer scanlines, 1px height, 2px wider + // Data type for input is Float32 (4 bytes) + // keep only three scanlines in memory at a time, make room for initial and final nodata + std::unique_ptr scanLine = qgis::make_unique( inputBlock->dataType(), scanLineWidth, 1 ); + // Note: output block is not 2px wider and it is an image + // Prepare context and queue + cl::Context ctx = QgsOpenClUtils::context(); + cl::CommandQueue queue( ctx ); + + // Cast to float (because double just crashes on some GPUs) + std::vector rasterParams; + rasterParams.push_back( inputBlock->noDataValue() ); + rasterParams.push_back( outputBlock->noDataValue() ); + rasterParams.push_back( mZFactor ); + rasterParams.push_back( cellXSize ); + rasterParams.push_back( cellYSize ); + rasterParams.push_back( static_cast( mOpacity ) ); // 5 + + // For fast formula from GDAL DEM + rasterParams.push_back( cos_az_mul_cos_alt_mul_z_mul_254 ); // 6 + rasterParams.push_back( sin_az_mul_cos_alt_mul_z_mul_254 ); // 7 + rasterParams.push_back( square_z ); // 8 + rasterParams.push_back( sin_altRadians_mul_254 ); // 9 + + // For multidirectional fast formula + rasterParams.push_back( sin_altRadians_mul_127 ); // 10 + rasterParams.push_back( cos225_az_mul_cos_alt_mul_z_mul_127 ); // 11 + rasterParams.push_back( cos_alt_mul_z_mul_127 ); // 12 + + // Default color for nodata (BGR components) + rasterParams.push_back( static_cast( qBlue( defaultNodataColor ) ) ); // 13 + rasterParams.push_back( static_cast( qGreen( defaultNodataColor ) ) ); // 14 + rasterParams.push_back( static_cast( qRed( defaultNodataColor ) ) ); // 15 + rasterParams.push_back( static_cast( qAlpha( defaultNodataColor ) ) / 255.0f ); // 16 + + // Whether use multidirectional + rasterParams.push_back( static_cast( mMultiDirectional ) ); // 17 + + + cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr ); + cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr ); + cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr ); + cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr ); + cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer}; + // Note that result buffer is an image + cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, outputDataTypeSize * width, nullptr, nullptr ); + + static cl::Program program; + static std::once_flag programBuilt; + std::call_once( programBuilt, [ = ]() { - outputBlock->setColor( i, j, myDefaultColor ); - continue; - } + // Create a program from the kernel source + program = QgsOpenClUtils::buildProgram( ctx, source, QgsOpenClUtils::ExceptionBehavior::Throw ); + } ); - qgssize iUp, iDown, jLeft, jRight; - if ( i == 0 ) + // Disable program cache when developing and testing cl program + // program = QgsOpenClUtils::buildProgram( ctx, source, QgsOpenClUtils::ExceptionBehavior::Throw ); + + // Create the OpenCL kernel + auto kernel = cl::KernelFunctor < + cl::Buffer &, + cl::Buffer &, + cl::Buffer &, + cl::Buffer &, + cl::Buffer & + > ( program, "processNineCellWindow" ); + + + // Rotating buffer index + std::vector rowIndex = {0, 1, 2}; + + for ( int i = 0; i < height; i++ ) { - iUp = i; - iDown = i + 1; - } - else if ( i < ( qgssize )height - 1 ) - { - iUp = i - 1; - iDown = i + 1; - } - else - { - iUp = i - 1; - iDown = i; - } - - if ( j == 0 ) - { - jLeft = j; - jRight = j + 1; - } - else if ( j < ( qgssize )width - 1 ) - { - jLeft = j - 1; - jRight = j + 1; - } - else - { - jLeft = j - 1; - jRight = j; - } - - float x11; - float x21; - float x31; - float x12; - float x22; // Working cell - float x32; - float x13; - float x23; - float x33; - - // This is center cell. It is not nodata. Use this in place of nodata neighbors - x22 = inputBlock->value( i, j ); - - x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft ); - x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft ); - x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft ); - - x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j ); - // x22 - x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j ); - - x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight ); - x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight ); - x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight ); - - float derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize ); - float derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize ); - - float grayValue; - if ( !mMultiDirectional ) - { - // Standard single direction hillshade - // Fast formula from GDAL DEM - grayValue = qBound( 0.0f, ( sinZenithRad_mul_254 - - ( derY * cos_az_mul_cos_alt_mul_z_mul_254 - - derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) / - std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) ) - , 255.0f ); - } - else - { - // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf - // Fast formula from GDAL DEM - const float xx = derX * derX; - const float yy = derY * derY; - const float xx_plus_yy = xx + yy; - // Flat? - if ( xx_plus_yy == 0.0 ) + if ( feedback && feedback->isCanceled() ) { - grayValue = qBound( 0.0f, static_cast( 1.0 + sinZenithRad_mul_254 ), 255.0f ); + break; + } + + if ( feedback ) + { + feedback->setProgress( 100.0 * static_cast< double >( i ) / height ); + } + + if ( i == 0 ) + { + // Fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row + scanLine->resetNoDataValue(); + queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); + // Read first row + memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i, 0 ), inputSize ); + queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); // row 0 + // Second row + memcpy( scanLine->bits( 0, 1 ), inputBlock->bits( i + 1, 0 ), inputSize ); + queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine->bits( ) ); // } else { - // ... then the shade value from different azimuth - float val225_mul_127 = sinZenithRad_mul_127 + - ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127; - val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127; - float val270_mul_127 = sinZenithRad_mul_127 - - derX * cos_alt_mul_z_mul_127; - val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127; - float val315_mul_127 = sinZenithRad_mul_127 + - ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127; - val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127; - float val360_mul_127 = sinZenithRad_mul_127 - - derY * cos_alt_mul_z_mul_127; - val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127; - - // ... then the weighted shading - const float weight_225 = 0.5 * xx_plus_yy - derX * derY; - const float weight_270 = xx; - const float weight_315 = xx_plus_yy - weight_225; - const float weight_360 = yy; - const float cang_mul_127 = ( - ( weight_225 * val225_mul_127 + - weight_270 * val270_mul_127 + - weight_315 * val315_mul_127 + - weight_360 * val360_mul_127 ) / xx_plus_yy ) / - ( 1 + square_z * xx_plus_yy ); - - grayValue = qBound( 0.0f, 1.0f + cang_mul_127, 255.0f ); + // Normally fetch only scanLine3 and move forward one row + // Read scanline 3, fill the last row with nodata values if it's the last iteration + if ( i == inputBlock->height() - 1 ) + { + scanLine->resetNoDataValue(); + queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine->bits( ) ); + } + else // Overwrite from input, skip first and last + { + queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, inputDataTypeSize * 1 /* offset 1 */, inputSize, inputBlock->bits( i + 1, 0 ) ); + } } - } - float currentAlpha = mOpacity; - if ( mRasterTransparency ) - { - currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0; - } - if ( mAlphaBand > 0 ) - { - currentAlpha *= alphaBlock->value( i ) / 255.0; - } + kernel( cl::EnqueueArgs( + queue, + cl::NDRange( width ) + ), + *scanLineBuffer[rowIndex[0]], + *scanLineBuffer[rowIndex[1]], + *scanLineBuffer[rowIndex[2]], + resultLineBuffer, + rasterParamsBuffer + ); - if ( qgsDoubleNear( currentAlpha, 1.0 ) ) - { - outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) ); - } - else - { - outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) ); + queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, outputDataTypeSize * outputBlock->width( ), outputBlock->bits( i, 0 ) ); + std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() ); } } - } + catch ( cl::Error &e ) + { + QgsMessageLog::logMessage( QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what( ) ).arg( QgsOpenClUtils::errorText( e.err( ) ) ), + QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical ); + QgsOpenClUtils::setEnabled( false ); + QgsMessageLog::logMessage( QObject::tr( "OpenCL has been disabled, you can re-enable it in the options dialog." ), + QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical ); + } + + } // End of OpenCL processing path + else // Use the CPU and the original algorithm + { + +#endif + + for ( qgssize i = 0; i < static_cast( height ); i++ ) + { + + for ( qgssize j = 0; j < static_cast( width ); j++ ) + { + + if ( inputBlock->isNoData( i, j ) ) + { + outputBlock->setColor( static_cast( i ), static_cast( j ), defaultNodataColor ); + continue; + } + + qgssize iUp, iDown, jLeft, jRight; + if ( i == 0 ) + { + iUp = i; + iDown = i + 1; + } + else if ( i < static_cast( height ) - 1 ) + { + iUp = i - 1; + iDown = i + 1; + } + else + { + iUp = i - 1; + iDown = i; + } + + if ( j == 0 ) + { + jLeft = j; + jRight = j + 1; + } + else if ( j < static_cast( width ) - 1 ) + { + jLeft = j - 1; + jRight = j + 1; + } + else + { + jLeft = j - 1; + jRight = j; + } + + double x11; + double x21; + double x31; + double x12; + double x22; // Working cell + double x32; + double x13; + double x23; + double x33; + + // This is center cell. It is not nodata. Use this in place of nodata neighbors + x22 = inputBlock->value( i, j ); + + x11 = inputBlock->isNoData( iUp, jLeft ) ? x22 : inputBlock->value( iUp, jLeft ); + x21 = inputBlock->isNoData( i, jLeft ) ? x22 : inputBlock->value( i, jLeft ); + x31 = inputBlock->isNoData( iDown, jLeft ) ? x22 : inputBlock->value( iDown, jLeft ); + + x12 = inputBlock->isNoData( iUp, j ) ? x22 : inputBlock->value( iUp, j ); + // x22 + x32 = inputBlock->isNoData( iDown, j ) ? x22 : inputBlock->value( iDown, j ); + + x13 = inputBlock->isNoData( iUp, jRight ) ? x22 : inputBlock->value( iUp, jRight ); + x23 = inputBlock->isNoData( i, jRight ) ? x22 : inputBlock->value( i, jRight ); + x33 = inputBlock->isNoData( iDown, jRight ) ? x22 : inputBlock->value( iDown, jRight ); + + double derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellXSize ); + double derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33, cellYSize ); + + // Fast formula + + double grayValue; + if ( !mMultiDirectional ) + { + // Standard single direction hillshade + grayValue = qBound( 0.0, ( sin_altRadians_mul_254 - + ( derY * cos_az_mul_cos_alt_mul_z_mul_254 - + derX * sin_az_mul_cos_alt_mul_z_mul_254 ) ) / + std::sqrt( 1 + square_z * ( derX * derX + derY * derY ) ) + , 255.0 ); + } + else + { + // Weighted multi direction as in http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf + // Fast formula from GDAL DEM + const float xx = derX * derX; + const float yy = derY * derY; + const float xx_plus_yy = xx + yy; + // Flat? + if ( xx_plus_yy == 0.0 ) + { + grayValue = qBound( 0.0f, static_cast( 1.0 + sin_altRadians_mul_254 ), 255.0f ); + } + else + { + // ... then the shade value from different azimuth + float val225_mul_127 = sin_altRadians_mul_127 + + ( derX - derY ) * cos225_az_mul_cos_alt_mul_z_mul_127; + val225_mul_127 = ( val225_mul_127 <= 0.0 ) ? 0.0 : val225_mul_127; + float val270_mul_127 = sin_altRadians_mul_127 - + derX * cos_alt_mul_z_mul_127; + val270_mul_127 = ( val270_mul_127 <= 0.0 ) ? 0.0 : val270_mul_127; + float val315_mul_127 = sin_altRadians_mul_127 + + ( derX + derY ) * cos225_az_mul_cos_alt_mul_z_mul_127; + val315_mul_127 = ( val315_mul_127 <= 0.0 ) ? 0.0 : val315_mul_127; + float val360_mul_127 = sin_altRadians_mul_127 - + derY * cos_alt_mul_z_mul_127; + val360_mul_127 = ( val360_mul_127 <= 0.0 ) ? 0.0 : val360_mul_127; + + // ... then the weighted shading + const float weight_225 = 0.5 * xx_plus_yy - derX * derY; + const float weight_270 = xx; + const float weight_315 = xx_plus_yy - weight_225; + const float weight_360 = yy; + const float cang_mul_127 = ( + ( weight_225 * val225_mul_127 + + weight_270 * val270_mul_127 + + weight_315 * val315_mul_127 + + weight_360 * val360_mul_127 ) / xx_plus_yy ) / + ( 1 + square_z * xx_plus_yy ); + + grayValue = qBound( 0.0f, 1.0f + cang_mul_127, 255.0f ); + } + } + + double currentAlpha = mOpacity; + if ( mRasterTransparency ) + { + currentAlpha = mRasterTransparency->alphaValue( x22, mOpacity * 255 ) / 255.0; + } + if ( mAlphaBand > 0 ) + { + currentAlpha *= alphaBlock->value( i ) / 255.0; + } + + if ( qgsDoubleNear( currentAlpha, 1.0 ) ) + { + outputBlock->setColor( i, j, qRgba( grayValue, grayValue, grayValue, 255 ) ); + } + else + { + outputBlock->setColor( i, j, qRgba( currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * grayValue, currentAlpha * 255 ) ); + } + } + } + +#ifdef HAVE_OPENCL + } // End of switch in case OpenCL is not available or enabled #ifdef QGISDEBUG if ( QgsSettings().value( QStringLiteral( "Map/logCanvasRefreshEvent" ), false ).toBool() ) { - QgsMessageLog::logMessage( QStringLiteral( "CPU processing time for hillshade (%2 x %3 ): %4 ms" ) + QgsMessageLog::logMessage( QStringLiteral( "%1 processing time for hillshade (%2 x %3 ): %4 ms" ) + .arg( useOpenCL ? QStringLiteral( "OpenCL" ) : QStringLiteral( "CPU" ) ) .arg( width ) .arg( height ) .arg( std::chrono::duration_cast( std::chrono::system_clock::now() - startTime ).count() ), tr( "Rendering" ) ); } +#endif #endif return outputBlock.release(); diff --git a/src/core/raster/qgsrasterblock.h b/src/core/raster/qgsrasterblock.h index 6b3fd8c2969..6e1b10b4486 100644 --- a/src/core/raster/qgsrasterblock.h +++ b/src/core/raster/qgsrasterblock.h @@ -238,16 +238,26 @@ class CORE_EXPORT QgsRasterBlock * \param row row index * \param column column index * \returns true if value is no data */ - bool isNoData( int row, int column ) + bool isNoData( int row, int column ) const { return isNoData( static_cast< qgssize >( row ) * mWidth + column ); } + /** + * \brief Check if value at position is no data + * \param row row index + * \param column column index + * \returns true if value is no data */ + bool isNoData( qgssize row, qgssize column ) const + { + return isNoData( row * static_cast< qgssize >( mWidth ) + column ); + } + /** * \brief Check if value at position is no data * \param index data matrix index (long type in Python) * \returns true if value is no data */ - bool isNoData( qgssize index ) + bool isNoData( qgssize index ) const { if ( !mHasNoDataValue && !mNoDataBitmap ) return false; diff --git a/src/ui/qgsoptionsbase.ui b/src/ui/qgsoptionsbase.ui index a138dd21673..6d9dd07ea4f 100644 --- a/src/ui/qgsoptionsbase.ui +++ b/src/ui/qgsoptionsbase.ui @@ -287,6 +287,18 @@ :/images/themes/default/mIconWarning.svg:/images/themes/default/mIconWarning.svg + + + Acceleration + + + Configure GPU for processing algorithms + + + + :/images/themes/default/mIconGPU.svg:/images/themes/default/mIconGPU.svg + + @@ -320,7 +332,7 @@ - 0 + 16 @@ -349,7 +361,7 @@ 0 0 - 843 + 556 885 @@ -1071,7 +1083,7 @@ 0 0 544 - 1094 + 1096 @@ -3204,7 +3216,7 @@ 0 0 613 - 636 + 646 @@ -5307,6 +5319,61 @@ The bigger the number, the faster zooming with the mouse wheel will be. + + + + + + <html><head/><body><p>Some of the internal C++ processing core algorithms and renderers can take advantage of an OpenCL compatible device to increase the performances.<br/><span style=" font-weight:600;">QGIS OpenCL support is highly experimental and can crash QGIS because of bugs in the underlying libraries, enable at your own risk!</span></p></body></html> + + + true + + + + + + + The following OpenCL devices were found on this system (changing the default device requires QGIS to be restarted). + + + + + + + + + + <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN" "http://www.w3.org/TR/REC-html40/strict.dtd"> +<html><head><meta name="qrichtext" content="1" /><style type="text/css"> +p, li { white-space: pre-wrap; } +</style></head><body style=" font-family:'Noto Sans'; font-size:9pt; font-weight:400; font-style:normal;"> +<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">Placemark for OpenCL information results (mGPUInfoTextBrowser)</p></body></html> + + + + + + + Enable OpenCL acceleration + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + @@ -5672,34 +5739,6 @@ The bigger the number, the faster zooming with the mouse wheel will be. - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/src/analysis/testqgsninecellfilters.cpp b/tests/src/analysis/testqgsninecellfilters.cpp index 89460c96640..b6d300ddc38 100644 --- a/tests/src/analysis/testqgsninecellfilters.cpp +++ b/tests/src/analysis/testqgsninecellfilters.cpp @@ -21,9 +21,16 @@ #include "qgsruggednessfilter.h" #include "qgstotalcurvaturefilter.h" #include "qgsapplication.h" +#include "qgssettings.h" + +#ifdef HAVE_OPENCL +#include "qgsopenclutils.h" +#endif #include +// If true regenerate raster reference images +const bool REGENERATE_REFERENCES = false; class TestNineCellFilters : public QObject { @@ -33,15 +40,25 @@ class TestNineCellFilters : public QObject private slots: void initTestCase(); + void init(); + + void testHillshade(); void testSlope(); void testAspect(); - void testHillshade(); void testRuggedness(); void testTotalCurvature(); +#ifdef HAVE_OPENCL + void testHillshadeCl(); + void testSlopeCl(); + void testAspectCl(); + void testRuggednessCl(); +#endif private: - template void _testAlg( const QString &name ); + void _rasterCompare( QgsAlignRaster::RasterInfo &out, QgsAlignRaster::RasterInfo &ref ); + + template void _testAlg( const QString &name, bool useOpenCl = false ); static QString referenceFile( const QString &name ) { @@ -55,6 +72,13 @@ class TestNineCellFilters : public QObject }; +void TestNineCellFilters::init() +{ +#ifdef HAVE_OPENCL + // Reset to default in case some tests mess it up + QgsOpenClUtils::setSourcePath( QDir( QgsApplication::pkgDataPath() ).absoluteFilePath( QStringLiteral( "resources/opencl_programs" ) ) ); +#endif +} void TestNineCellFilters::initTestCase() { @@ -65,36 +89,37 @@ void TestNineCellFilters::initTestCase() } template -void TestNineCellFilters::_testAlg( const QString &name ) +void TestNineCellFilters::_testAlg( const QString &name, bool useOpenCl ) { +#ifdef HAVE_OPENCL + QgsOpenClUtils::setEnabled( useOpenCl ); + QString tmpFile( tempFile( name + ( useOpenCl ? "_opencl" : "" ) ) ); +#else QString tmpFile( tempFile( name ) ); +#endif QString refFile( referenceFile( name ) ); - - QgsAlignRaster::RasterInfo in( SRC_FILE ); - QSize inSize( in.rasterSize() ); - QSizeF inCellSize( in.cellSize( ) ); - T slopefilter( SRC_FILE, tmpFile, "GTiff" ); - int res = slopefilter.processRaster(); + T ninecellFilter( SRC_FILE, tmpFile, "GTiff" ); + int res = ninecellFilter.processRaster(); QVERIFY( res == 0 ); // Produced file QgsAlignRaster::RasterInfo out( tmpFile ); QVERIFY( out.isValid() ); - // Reference file + // Regenerate reference rasters + if ( ! useOpenCl && REGENERATE_REFERENCES ) + { + if ( QFile::exists( refFile ) ) + { + QFile::remove( refFile ); + } + QVERIFY( QFile::copy( tmpFile, refFile ) ); + } + + // Reference QgsAlignRaster::RasterInfo ref( refFile ); - QSize refSize( ref.rasterSize() ); - QSizeF refCellSize( ref.cellSize( ) ); - - QCOMPARE( out.rasterSize(), inSize ); - QCOMPARE( out.cellSize(), inCellSize ); - QCOMPARE( out.rasterSize(), refSize ); - QCOMPARE( out.cellSize(), refCellSize ); - - double refId1( ref.identify( 4081812, 2431750 ) ); - double refId2( ref.identify( 4081312, 2431350 ) ); - QCOMPARE( out.identify( 4081812, 2431750 ), refId1 ); - QCOMPARE( out.identify( 4081312, 2431350 ), refId2 ); + //qDebug() << "Comparing " << tmpFile << refFile; + _rasterCompare( out, ref ); } @@ -104,30 +129,93 @@ void TestNineCellFilters::testSlope() _testAlg( QStringLiteral( "slope" ) ); } - void TestNineCellFilters::testAspect() { _testAlg( QStringLiteral( "aspect" ) ); } +#ifdef HAVE_OPENCL +void TestNineCellFilters::testSlopeCl() +{ + _testAlg( QStringLiteral( "slope" ), true ); +} + +void TestNineCellFilters::testAspectCl() +{ + _testAlg( QStringLiteral( "aspect" ), true ); +} + +void TestNineCellFilters::testHillshadeCl() +{ + _testAlg( QStringLiteral( "hillshade" ), true ); +} + +void TestNineCellFilters::testRuggednessCl() +{ + _testAlg( QStringLiteral( "ruggedness" ), true ); +} + +#endif void TestNineCellFilters::testHillshade() { _testAlg( QStringLiteral( "hillshade" ) ); } - void TestNineCellFilters::testRuggedness() { _testAlg( QStringLiteral( "ruggedness" ) ); } +void TestNineCellFilters::_rasterCompare( QgsAlignRaster::RasterInfo &out, QgsAlignRaster::RasterInfo &ref ) +{ + QSize refSize( ref.rasterSize() ); + QSizeF refCellSize( ref.cellSize( ) ); + QgsAlignRaster::RasterInfo in( SRC_FILE ); + QSize inSize( in.rasterSize() ); + QSizeF inCellSize( in.cellSize( ) ); + QCOMPARE( out.rasterSize(), inSize ); + QCOMPARE( out.cellSize(), inCellSize ); + QCOMPARE( out.rasterSize(), refSize ); + QCOMPARE( out.cellSize(), refCellSize ); + + // If the values differ less than tolerance they are considered equal + double tolerance = 0.0001; + + // Check three points + std::map controlPoints; + controlPoints[4081812] = 2431750; + controlPoints[4081312] = 2431350; + controlPoints[4080263] = 2429558; + // South West corner + controlPoints[4081512] = 2431550; + // North east corner + controlPoints[4085367] = 2434940; + // North west corner + controlPoints[4078263] = 2434936; + // South east corner + controlPoints[4085374] = 2428551; + + for ( const auto &cp : controlPoints ) + { + int x = cp.first; + int y = cp.second; + double outVal = out.identify( x, y ); + double refVal = ref.identify( x, y ); + double diff( qAbs( outVal - refVal ) ); + //qDebug() << outVal << refVal; + //qDebug() << "Identify " << x << "," << y << " diff " << diff << " check: < " << tolerance; + QVERIFY( diff <= tolerance ); + } + +} void TestNineCellFilters::testTotalCurvature() { _testAlg( QStringLiteral( "totalcurvature" ) ); } + QGSTEST_MAIN( TestNineCellFilters ) #include "testqgsninecellfilters.moc" diff --git a/tests/src/core/CMakeLists.txt b/tests/src/core/CMakeLists.txt index fa4932b752d..9e040c82c0f 100755 --- a/tests/src/core/CMakeLists.txt +++ b/tests/src/core/CMakeLists.txt @@ -205,6 +205,14 @@ IF(WITH_QTWEBKIT) ) ENDIF(WITH_QTWEBKIT) + +IF(HAVE_OPENCL) + SET(TESTS ${TESTS} + testqgsopenclutils.cpp + ) +ENDIF(HAVE_OPENCL) + + FOREACH(TESTSRC ${TESTS}) ADD_QGIS_TEST(${TESTSRC}) ENDFOREACH(TESTSRC) diff --git a/tests/src/core/testqgsopenclutils.cpp b/tests/src/core/testqgsopenclutils.cpp new file mode 100644 index 00000000000..028959b1315 --- /dev/null +++ b/tests/src/core/testqgsopenclutils.cpp @@ -0,0 +1,242 @@ +/*************************************************************************** + testqgsopenclutils.cpp - TestQgsOpenClUtils + + --------------------- + begin : 11.4.2018 + copyright : (C) 2018 by Alessandro Pasotti + email : elpaso at itopen dot it + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include "qgstest.h" +#include +#include +#include +#include +#include + +//header for class being tested +#include +#include +#include + +class TestQgsOpenClUtils: public QObject +{ + Q_OBJECT + public: + + //void testRunMakeProgram(); + + private slots: + void initTestCase();// will be called before the first testfunction is executed. + void cleanupTestCase();// will be called after the last testfunction was executed. + void init();// will be called before each testfunction is executed. + void cleanup();// will be called after every testfunction. + + void TestEnable(); + void TestDisable(); + void TestAvailable(); + void testMakeRunProgram(); + void testProgramSource(); + void testContext(); + void testDevices(); + + // For benchmarking performance testing + void testHillshadeCPU(); + void testHillshadeGPU(); + + private: + + void _testMakeRunProgram(); + void _testMakeHillshade( const int loops ); + + cl::Program buildProgram( const cl::Context &context, const QString &source ) + { + cl::Program program( context, source.toStdString( ) ); + program.build( "-cl-std=CL1.1" ); + return program; + } + + std::string source() + { + std::string pgm = R"CL( + __kernel void vectorAdd(__global float *a, __global float *b, __global float *c) + { + const int id = get_global_id(0); + c[id] = a[id] + b[id]; + } + )CL"; + return pgm; + + } + + QgsRasterLayer *mFloat32RasterLayer = nullptr; +}; + + +void TestQgsOpenClUtils::init() +{ + // Reset to default in case some tests mess it up + QgsOpenClUtils::setSourcePath( QDir( QgsApplication::pkgDataPath() ).absoluteFilePath( QStringLiteral( "resources/opencl_programs" ) ) ); +} + +void TestQgsOpenClUtils::cleanup() +{ +} + +void TestQgsOpenClUtils::initTestCase() +{ + // Runs once before any tests are run + + // Set up the QgsSettings environment + QCoreApplication::setOrganizationName( QStringLiteral( "QGIS" ) ); + QCoreApplication::setOrganizationDomain( QStringLiteral( "qgis.org" ) ); + QCoreApplication::setApplicationName( QStringLiteral( "QGIS-TEST" ) ); + + QgsApplication::init(); + QgsApplication::initQgis(); + + QString float32FileName = QStringLiteral( TEST_DATA_DIR ) + '/' + "/raster/band1_float32_noct_epsg4326.tif"; + QFileInfo float32RasterFileInfo( float32FileName ); + mFloat32RasterLayer = new QgsRasterLayer( float32RasterFileInfo.filePath(), + float32RasterFileInfo.completeBaseName() ); +} + + +void TestQgsOpenClUtils::cleanupTestCase() +{ + // Runs once after all tests are run + QgsApplication::exitQgis(); +} + + +void TestQgsOpenClUtils::TestEnable() +{ + QgsOpenClUtils::setEnabled( true ); + QVERIFY( QgsOpenClUtils::enabled() ); +} + +void TestQgsOpenClUtils::TestDisable() +{ + QgsOpenClUtils::setEnabled( false ); + QVERIFY( !QgsOpenClUtils::enabled() ); +} + +void TestQgsOpenClUtils::TestAvailable() +{ + QVERIFY( QgsOpenClUtils::available() ); +} + + +void TestQgsOpenClUtils::testMakeRunProgram() +{ + // Run twice to check for valid command queue in the second call + _testMakeRunProgram(); + _testMakeRunProgram(); +} + +void TestQgsOpenClUtils::_testMakeRunProgram() +{ + + cl_int err = 0; + + QVERIFY( err == 0 ); + + cl::Context ctx = QgsOpenClUtils::context(); + cl::CommandQueue queue( ctx ); + + std::vector a_vec = {1, 10, 100}; + std::vector b_vec = {1, 10, 100}; + std::vector c_vec = {-1, -1, -1}; + cl::Buffer a_buf( queue, a_vec.begin(), a_vec.end(), true ); + cl::Buffer b_buf( queue, b_vec.begin(), b_vec.end(), true ); + cl::Buffer c_buf( queue, c_vec.begin(), c_vec.end(), false ); + + cl::Program program = QgsOpenClUtils::buildProgram( ctx, QString::fromStdString( source() ) ); + + auto kernel = + cl::KernelFunctor < + cl::Buffer &, + cl::Buffer &, + cl::Buffer & + > ( program, "vectorAdd" ); + + kernel( cl::EnqueueArgs( + queue, + cl::NDRange( 3 ) + ), + a_buf, + b_buf, + c_buf + ); + + cl::copy( queue, c_buf, c_vec.begin(), c_vec.end() ); + for ( size_t i = 0; i < c_vec.size(); ++i ) + { + QCOMPARE( c_vec[i], a_vec[i] + b_vec[i] ); + } +} + +void TestQgsOpenClUtils::testProgramSource() +{ + QgsOpenClUtils::setSourcePath( QDir::tempPath() ); + QTemporaryFile tmpFile( QDir::tempPath() + "/XXXXXX.cl" ); + tmpFile.open( ); + tmpFile.write( QByteArray::fromStdString( source( ) ) ); + tmpFile.flush(); + QString baseName = tmpFile.fileName().replace( ".cl", "" ).replace( QDir::tempPath(), "" ); + QVERIFY( ! QgsOpenClUtils::sourceFromBaseName( baseName ).isEmpty() ); +} + +void TestQgsOpenClUtils::testContext() +{ + QVERIFY( QgsOpenClUtils::context()() != nullptr ); +} + +void TestQgsOpenClUtils::testDevices() +{ + std::vector _devices( QgsOpenClUtils::devices( ) ); + QVERIFY( _devices.size() > 0 ); + qDebug() << QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Name, _devices.at( 0 ) ); + qDebug() << QgsOpenClUtils::deviceInfo( QgsOpenClUtils::Info::Type, _devices.at( 0 ) ); +} + +void TestQgsOpenClUtils::_testMakeHillshade( const int loops ) +{ + for ( int i = 0 ; i < loops; i++ ) + { + QgsHillshadeRenderer renderer( mFloat32RasterLayer->dataProvider(), 1, 35.0, 5000.0 ); + // Note: CPU time grows linearly with raster dimensions while GPU time is roughly constant + // 200x200 px gives even times on my testing machine + renderer.block( 0, mFloat32RasterLayer->extent(), 200, 200 ); + } +} + +void TestQgsOpenClUtils::testHillshadeGPU() +{ + QgsOpenClUtils::setEnabled( true ); + QBENCHMARK + { + _testMakeHillshade( 1 ); + } +} + +void TestQgsOpenClUtils::testHillshadeCPU() +{ + QgsOpenClUtils::setEnabled( false ); + QBENCHMARK + { + _testMakeHillshade( 1 ); + } +} + + + +QGSTEST_MAIN( TestQgsOpenClUtils ) +#include "testqgsopenclutils.moc" diff --git a/tests/testdata/analysis/hillshade.tif b/tests/testdata/analysis/hillshade.tif index 92d0832d3b3..96b1154b5af 100644 Binary files a/tests/testdata/analysis/hillshade.tif and b/tests/testdata/analysis/hillshade.tif differ diff --git a/tests/testdata/opencl/aspect.tif b/tests/testdata/opencl/aspect.tif new file mode 100644 index 00000000000..b901524ac19 Binary files /dev/null and b/tests/testdata/opencl/aspect.tif differ diff --git a/tests/testdata/opencl/dem.tif b/tests/testdata/opencl/dem.tif new file mode 100644 index 00000000000..8bf7cd5d01b Binary files /dev/null and b/tests/testdata/opencl/dem.tif differ diff --git a/tests/testdata/opencl/slope.tif b/tests/testdata/opencl/slope.tif new file mode 100644 index 00000000000..9c0a36f3658 Binary files /dev/null and b/tests/testdata/opencl/slope.tif differ