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

[feature] OpenCL support
This commit is contained in:
Alessandro Pasotti 2018-08-08 16:33:48 +02:00 committed by GitHub
commit 0b502ff5b9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
45 changed files with 2780 additions and 359 deletions

View File

@ -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)

View File

@ -57,6 +57,8 @@
#cmakedefine HAVE_SERVER_PYTHON_PLUGINS
#cmakedefine HAVE_OPENCL
#cmakedefine ENABLE_MODELTEST
#cmakedefine HAVE_3D

View File

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

View File

@ -0,0 +1,14 @@
<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 22 22">
<defs id="defs3051">
<style type="text/css" id="current-color-scheme">
.ColorScheme-Text {
color:#4d4d4d;
}
</style>
</defs>
<path
style="fill:currentColor;fill-opacity:1;stroke:none"
d="M 6 3 L 6 5 L 5 5 L 5 6 L 3 6 L 3 7 L 5 7 L 5 9 L 3 9 L 3 10 L 5 10 L 5 12 L 3 12 L 3 13 L 5 13 L 5 15 L 3 15 L 3 16 L 5 16 L 5 17 L 6 17 L 6 19 L 7 19 L 7 17 L 9 17 L 9 19 L 10 19 L 10 17 L 12 17 L 12 19 L 13 19 L 13 17 L 15 17 L 15 19 L 16 19 L 16 17 L 17 17 L 17 16 L 19 16 L 19 15 L 17 15 L 17 13 L 19 13 L 19 12 L 17 12 L 17 10 L 19 10 L 19 9 L 17 9 L 17 7 L 19 7 L 19 6 L 17 6 L 17 5 L 16 5 L 16 3 L 15 3 L 15 5 L 13 5 L 13 3 L 12 3 L 12 5 L 10 5 L 10 3 L 9 3 L 9 5 L 7 5 L 7 3 L 6 3 z "
class="ColorScheme-Text"
/>
</svg>

After

Width:  |  Height:  |  Size: 832 B

View File

@ -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
};
/************************************************************************

View File

@ -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:

View File

@ -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
};
/************************************************************************

View File

@ -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

View File

@ -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)

View File

@ -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;
}
}
}

View File

@ -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);
}

View File

@ -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;
}

View File

@ -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 ) ) ) );
}
}
}

View File

@ -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]);
}
}
}

View File

@ -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);
}
}

View File

@ -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;
}
}
}

View File

@ -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);
}

View File

@ -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(

View File

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

View File

@ -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

View File

@ -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<float>( lightAzimuth ) )
, mLightAngle( static_cast<float>( lightAngle ) )
, mCosZenithRad( std::cos( static_cast<float>( lightAngle * M_PI ) / 180.0f ) )
, mSinZenithRad( std::sin( static_cast<float>( lightAngle * M_PI ) / 180.0f ) )
, mAzimuthRad( static_cast<float>( lightAzimuth * M_PI ) / 180.0f )
{
}
@ -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<float>( M_PI ) / 180.0f;
}
void QgsHillshadeFilter::setLightAngle( float angle )
{
mLightAngle = angle;
mCosZenithRad = std::cos( angle * static_cast<float>( M_PI ) / 180.0f );
mSinZenithRad = std::sin( angle * static_cast<float>( M_PI ) / 180.0f );
}
#ifdef HAVE_OPENCL
void QgsHillshadeFilter::addExtraRasterParams( std::vector<float> &params )
{
params.push_back( mCosZenithRad ); // cos_zenith_rad 5
params.push_back( mSinZenithRad ); // sin_zenith_rad 6
params.push_back( mAzimuthRad ); // azimuth_rad 7
}
#endif

View File

@ -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<float> &params ) override;
#endif
};
#endif // QGSHILLSHADEFILTER_H

View File

@ -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 <QFile>
#include <QDebug>
#include <QFileInfo>
#include <iterator>
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<float> scanLine( xSize + 2 );
QgsOpenClUtils::CPLAllocator<float> resultLine( xSize );
// Cast to float (because double just crashes on some GPUs)
std::vector<float> 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<int> 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;
}

View File

@ -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<float> &params )
{
Q_UNUSED( params );
}
virtual const QString openClProgramBaseName() const
{
return QString();
}
#endif
protected:
QString mInputFile;

View File

@ -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 )

View File

@ -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

View File

@ -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 );

View File

@ -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

View File

@ -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...
/////////////////////////////////////////////////////////////////////

View File

@ -97,6 +97,10 @@
#include "qgsgui.h"
#include "qgsnative.h"
#ifdef HAVE_OPENCL
#include "qgsopenclutils.h"
#endif
#include <QNetworkReply>
#include <QNetworkProxy>
#include <QAuthenticator>
@ -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() );

View File

@ -51,6 +51,10 @@
#include "qgslocatorwidget.h"
#include "qgslocatoroptionswidget.h"
#ifdef HAVE_OPENCL
#include "qgsopenclutils.h"
#endif
#include <QInputDialog>
#include <QFileDialog>
#include <QColorDialog>
@ -1076,6 +1080,44 @@ QgsOptions::QgsOptions( QWidget *parent, Qt::WindowFlags fl, const QList<QgsOpti
mOptionsStackedWidget->addWidget( 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<void( int )> 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.<br>"
"You may need to install additional libraries in order to enable OpenCL.<br>"
"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();

View File

@ -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)

530
src/core/qgsopenclutils.cpp Normal file
View File

@ -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 <QTextStream>
#include <QFile>
#include <QDebug>
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<cl::Device> QgsOpenClUtils::devices()
{
std::vector<cl::Platform> platforms;
cl::Platform::get( &platforms );
std::vector<cl::Device> existingDevices;
for ( auto &p : platforms )
{
std::string platver = p.getInfo<CL_PLATFORM_VERSION>();
QgsDebugMsg( QStringLiteral( "Found OpenCL platform %1: %2" ).arg( QString::fromStdString( platver ), QString::fromStdString( p.getInfo<CL_PLATFORM_NAME>() ) ) );
if ( platver.find( "OpenCL 1." ) != std::string::npos )
{
std::vector<cl::Device> _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<CL_PLATFORM_NAME>() ) ), 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<CL_DEVICE_VENDOR>() );
case Info::Profile:
return QString::fromStdString( device.getInfo<CL_DEVICE_PROFILE>() );
case Info::Version:
return QString::fromStdString( device.getInfo<CL_DEVICE_VERSION>() );
case Info::ImageSupport:
return device.getInfo<CL_DEVICE_IMAGE_SUPPORT>() ? QStringLiteral( "True" ) : QStringLiteral( "False" );
case Info::Image2dMaxHeight:
return QString::number( device.getInfo<CL_DEVICE_IMAGE2D_MAX_HEIGHT>() );
case Info::MaxMemAllocSize:
return QString::number( device.getInfo<CL_DEVICE_MAX_MEM_ALLOC_SIZE>() );
case Info::Image2dMaxWidth:
return QString::number( device.getInfo<CL_DEVICE_IMAGE2D_MAX_WIDTH>() );
case Info::Type:
{
unsigned long type( device.getInfo<CL_DEVICE_TYPE>() );
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<QgsOpenClUtils::HardwareType>();
return metaEnum.valueToKey( mappedType );
}
case Info::Name:
return QString::fromStdString( device.getInfo<CL_DEVICE_NAME>() );
}
}
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<cl::Platform> 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<CL_PLATFORM_VERSION>();
QgsDebugMsg( QStringLiteral( "Found OpenCL platform %1: %2" ).arg( QString::fromStdString( platver ), QString::fromStdString( p.getInfo<CL_PLATFORM_NAME>() ) ) );
if ( platver.find( "OpenCL 1." ) != std::string::npos )
{
std::vector<cl::Device> 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>() == 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>() == 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<CL_PLATFORM_NAME>() ) ) );
}
}
}
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<CL_DEVICE_NAME>() ) ),
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: <b>%9</b><br>"
"Name: <b>%1</b><br>"
"Vendor: <b>%2</b><br>"
"Profile: <b>%3</b><br>"
"Version: <b>%4</b><br>"
"Image support: <b>%5</b><br>"
"Max image2d width: <b>%6</b><br>"
"Max image2d height: <b>%7</b><br>"
"Max mem alloc size: <b>%8</b><br>"
).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;
}

289
src/core/qgsopenclutils.h Normal file
View File

@ -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 <CL/cl2.hpp>
#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<cl::Device> 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 <typename T>
struct CPLAllocator
{
public:
explicit CPLAllocator( unsigned long size ): mMem( static_cast<T *>( CPLMalloc( sizeof( T ) * size ) ) ) { }
~CPLAllocator()
{
CPLFree( static_cast<void *>( mMem ) );
}
void reset( T *newData )
{
if ( mMem )
CPLFree( static_cast<void *>( mMem ) );
mMem = newData;
}
void reset( unsigned long size )
{
reset( static_cast<T *>( 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

View File

@ -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

View File

@ -22,12 +22,15 @@
#include "qgsrasterinterface.h"
#include "qgsrasterblock.h"
#include "qgsrectangle.h"
#include "qgssettings.h"
#include "qgsmessagelog.h"
#include <memory>
#ifdef HAVE_OPENCL
#ifdef QGISDEBUG
#include "qgsmessagelog.h"
#include <chrono>
#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<float>( width );
float cellYSize = extent.height() / static_cast<float>( 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<float>( extent.width() ) / width;
float cellYSize = static_cast<float>( extent.height() ) / height;
float zenithRad = static_cast<float>( std::max( 0.0, 90 - mLightAngle ) * M_PI / 180.0 );
float azimuthRad = static_cast<float>( -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<float>( 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<float>( 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<std::chrono::system_clock> 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<QgsRasterBlock> scanLine = qgis::make_unique<QgsRasterBlock>( 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<float> 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<float>( 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<float>( qBlue( defaultNodataColor ) ) ); // 13
rasterParams.push_back( static_cast<float>( qGreen( defaultNodataColor ) ) ); // 14
rasterParams.push_back( static_cast<float>( qRed( defaultNodataColor ) ) ); // 15
rasterParams.push_back( static_cast<float>( qAlpha( defaultNodataColor ) ) / 255.0f ); // 16
// Whether use multidirectional
rasterParams.push_back( static_cast<float>( 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<int> 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<float>( 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<qgssize>( height ); i++ )
{
for ( qgssize j = 0; j < static_cast<qgssize>( width ); j++ )
{
if ( inputBlock->isNoData( i, j ) )
{
outputBlock->setColor( static_cast<int>( i ), static_cast<int>( j ), defaultNodataColor );
continue;
}
qgssize iUp, iDown, jLeft, jRight;
if ( i == 0 )
{
iUp = i;
iDown = i + 1;
}
else if ( i < static_cast<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 < static_cast<qgssize>( 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<float>( 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::milliseconds>( std::chrono::system_clock::now() - startTime ).count() ),
tr( "Rendering" ) );
}
#endif
#endif
return outputBlock.release();

View File

@ -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;

View File

@ -287,6 +287,18 @@
<normaloff>:/images/themes/default/mIconWarning.svg</normaloff>:/images/themes/default/mIconWarning.svg</iconset>
</property>
</item>
<item>
<property name="text">
<string>Acceleration</string>
</property>
<property name="toolTip">
<string>Configure GPU for processing algorithms</string>
</property>
<property name="icon">
<iconset resource="../../images/images.qrc">
<normaloff>:/images/themes/default/mIconGPU.svg</normaloff>:/images/themes/default/mIconGPU.svg</iconset>
</property>
</item>
</widget>
</item>
</layout>
@ -320,7 +332,7 @@
<item>
<widget class="QStackedWidget" name="mOptionsStackedWidget">
<property name="currentIndex">
<number>0</number>
<number>16</number>
</property>
<widget class="QWidget" name="mOptionsPageGeneral">
<layout class="QVBoxLayout" name="verticalLayout_3">
@ -349,7 +361,7 @@
<rect>
<x>0</x>
<y>0</y>
<width>843</width>
<width>556</width>
<height>885</height>
</rect>
</property>
@ -1071,7 +1083,7 @@
<x>0</x>
<y>0</y>
<width>544</width>
<height>1094</height>
<height>1096</height>
</rect>
</property>
<layout class="QVBoxLayout" name="verticalLayout_22">
@ -3204,7 +3216,7 @@
<x>0</x>
<y>0</y>
<width>613</width>
<height>636</height>
<height>646</height>
</rect>
</property>
<layout class="QVBoxLayout" name="verticalLayout_30">
@ -5307,6 +5319,61 @@ The bigger the number, the faster zooming with the mouse wheel will be.</string>
</item>
</layout>
</widget>
<widget class="QWidget" name="mOptionsPageAcceleration">
<layout class="QVBoxLayout" name="verticalLayout_29">
<item>
<widget class="QLabel" name="label_53">
<property name="text">
<string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Some of the internal C++ processing core algorithms and renderers can take advantage of an OpenCL compatible device to increase the performances.&lt;br/&gt;&lt;span style=&quot; font-weight:600;&quot;&gt;QGIS OpenCL support is highly experimental and can crash QGIS because of bugs in the underlying libraries, enable at your own risk!&lt;/span&gt;&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
</property>
<property name="wordWrap">
<bool>true</bool>
</property>
</widget>
</item>
<item>
<widget class="QLabel" name="label_64">
<property name="text">
<string>The following OpenCL devices were found on this system (changing the default device requires QGIS to be restarted).</string>
</property>
</widget>
</item>
<item>
<widget class="QComboBox" name="mOpenClDevicesCombo"/>
</item>
<item>
<widget class="QTextBrowser" name="mGPUInfoTextBrowser">
<property name="html">
<string>&lt;!DOCTYPE HTML PUBLIC &quot;-//W3C//DTD HTML 4.0//EN&quot; &quot;http://www.w3.org/TR/REC-html40/strict.dtd&quot;&gt;
&lt;html&gt;&lt;head&gt;&lt;meta name=&quot;qrichtext&quot; content=&quot;1&quot; /&gt;&lt;style type=&quot;text/css&quot;&gt;
p, li { white-space: pre-wrap; }
&lt;/style&gt;&lt;/head&gt;&lt;body style=&quot; font-family:'Noto Sans'; font-size:9pt; font-weight:400; font-style:normal;&quot;&gt;
&lt;p style=&quot; margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;&quot;&gt;Placemark for OpenCL information results (mGPUInfoTextBrowser)&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="mGPUEnableCheckBox">
<property name="text">
<string>Enable OpenCL acceleration</string>
</property>
</widget>
</item>
<item>
<spacer name="verticalSpacer_2">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
</layout>
</widget>
</widget>
</item>
</layout>
@ -5672,34 +5739,6 @@ The bigger the number, the faster zooming with the mouse wheel will be.</string>
</tabstops>
<resources>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
<include location="../../images/images.qrc"/>
</resources>
<connections>
<connection>

View File

@ -21,9 +21,16 @@
#include "qgsruggednessfilter.h"
#include "qgstotalcurvaturefilter.h"
#include "qgsapplication.h"
#include "qgssettings.h"
#ifdef HAVE_OPENCL
#include "qgsopenclutils.h"
#endif
#include <QDir>
// 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 <class T> void _testAlg( const QString &name );
void _rasterCompare( QgsAlignRaster::RasterInfo &out, QgsAlignRaster::RasterInfo &ref );
template <class T> 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 <class T>
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<QgsSlopeFilter>( QStringLiteral( "slope" ) );
}
void TestNineCellFilters::testAspect()
{
_testAlg<QgsAspectFilter>( QStringLiteral( "aspect" ) );
}
#ifdef HAVE_OPENCL
void TestNineCellFilters::testSlopeCl()
{
_testAlg<QgsSlopeFilter>( QStringLiteral( "slope" ), true );
}
void TestNineCellFilters::testAspectCl()
{
_testAlg<QgsAspectFilter>( QStringLiteral( "aspect" ), true );
}
void TestNineCellFilters::testHillshadeCl()
{
_testAlg<QgsHillshadeFilter>( QStringLiteral( "hillshade" ), true );
}
void TestNineCellFilters::testRuggednessCl()
{
_testAlg<QgsRuggednessFilter>( QStringLiteral( "ruggedness" ), true );
}
#endif
void TestNineCellFilters::testHillshade()
{
_testAlg<QgsHillshadeFilter>( QStringLiteral( "hillshade" ) );
}
void TestNineCellFilters::testRuggedness()
{
_testAlg<QgsRuggednessFilter>( 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<int, int> 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<QgsTotalCurvatureFilter>( QStringLiteral( "totalcurvature" ) );
}
QGSTEST_MAIN( TestNineCellFilters )
#include "testqgsninecellfilters.moc"

View File

@ -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)

View File

@ -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 <QObject>
#include <QString>
#include <QTemporaryFile>
#include <qgsapplication.h>
#include <chrono>
//header for class being tested
#include <qgsopenclutils.h>
#include <qgshillshaderenderer.h>
#include <qgsrasterlayer.h>
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<float> a_vec = {1, 10, 100};
std::vector<float> b_vec = {1, 10, 100};
std::vector<float> 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<cl::Device> _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"

Binary file not shown.

BIN
tests/testdata/opencl/aspect.tif vendored Normal file

Binary file not shown.

BIN
tests/testdata/opencl/dem.tif vendored Normal file

Binary file not shown.

BIN
tests/testdata/opencl/slope.tif vendored Normal file

Binary file not shown.