Skip to content

Commit dc4b1da

Browse files
elpasoelpaso
committedAug 8, 2018
[opencl] Fix small OpenCL alg issues
From comparison tests with CPU results + some minor speed improvements
·
ltr-3_40final-3_4_0
1 parent 30a62e1 commit dc4b1da

File tree

12 files changed

+111
-31
lines changed

12 files changed

+111
-31
lines changed
 

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

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,8 +30,6 @@ nodata value if not present or outside of the border. Must be implemented by sub
3030
%End
3131

3232

33-
34-
3533
};
3634

3735
/************************************************************************

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ Calculates output value from nine input values. The input values and the output
2929
nodata value if not present or outside of the border. Must be implemented by subclasses*
3030
%End
3131

32+
3233
};
3334

3435
/************************************************************************

‎resources/opencl_programs/hillshade.cl

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,11 @@
11
#include "calcfirstder.cl"
2+
<<<<<<< 30a62e1addd2541bdbecda4398381772eb8034d4
23

34
#ifndef M_PI
45
#define M_PI 3.14159265358979323846
56
#endif
7+
=======
8+
>>>>>>> [opencl] Fix small OpenCL alg issues
69

710
__kernel void processNineCellWindow( __global float *scanLine1,
811
__global float *scanLine2,
@@ -39,6 +42,7 @@ __kernel void processNineCellWindow( __global float *scanLine1,
3942
else
4043
{
4144

45+
<<<<<<< 30a62e1addd2541bdbecda4398381772eb8034d4
4246
float slope_rad = sqrt( derX * derX + derY * derY );
4347
slope_rad = atan( slope_rad );
4448
float aspect_rad;
@@ -50,6 +54,23 @@ __kernel void processNineCellWindow( __global float *scanLine1,
5054
resultLine[i] = max(0.0f, 255.0f * ( ( rasterParams[5] * cos( slope_rad ) ) +
5155
( rasterParams[6] * sin( slope_rad ) *
5256
cos( rasterParams[7] - aspect_rad ) ) ) );
57+
=======
58+
float zenith_rad = rasterParams[5];
59+
float azimuth_rad = rasterParams[6];
60+
float slope_rad = sqrt( derX * derX + derY * derY );
61+
62+
63+
slope_rad = atan( slope_rad );
64+
float aspect_rad;
65+
if ( derX == 0.0f && derY == 0.0f)
66+
aspect_rad = azimuth_rad / 2.0f;
67+
else
68+
aspect_rad = M_PI + atan2( derX, derY );
69+
70+
resultLine[i] = max(0.0f, 255.0f * ( ( cos( zenith_rad ) * cos( slope_rad ) ) +
71+
( sin( zenith_rad ) * sin( slope_rad ) *
72+
cos( azimuth_rad - aspect_rad ) ) ) );
73+
>>>>>>> [opencl] Fix small OpenCL alg issues
5374

5475
}
5576
}

‎src/analysis/raster/qgsalignraster.cpp

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -615,3 +615,4 @@ double QgsAlignRaster::RasterInfo::identify( double mx, double my )
615615

616616
return value;
617617
}
618+

‎src/analysis/raster/qgsaspectfilter.h

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -37,14 +37,15 @@ class ANALYSIS_EXPORT QgsAspectFilter: public QgsDerivativeFilter
3737
float *x13, float *x23, float *x33 ) override;
3838

3939

40-
41-
42-
// QgsNineCellFilter interface
40+
#ifdef HAVE_OPENCL
4341
private:
42+
4443
const QString openClProgramBaseName() const override
4544
{
4645
return QStringLiteral( "aspect" );
4746
}
47+
#endif
48+
4849
};
4950

5051
#endif // QGSASPECTFILTER_H

‎src/analysis/raster/qgshillshadefilter.cpp

Lines changed: 36 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -21,15 +21,19 @@
2121
QgsHillshadeFilter::QgsHillshadeFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat, double lightAzimuth,
2222
double lightAngle )
2323
: QgsDerivativeFilter( inputFile, outputFile, outputFormat )
24-
, mLightAzimuth( lightAzimuth )
25-
, mLightAngle( lightAngle )
24+
, mLightAzimuth( static_cast<float>( lightAzimuth ) )
25+
, mLightAngle( static_cast<float>( lightAngle ) )
26+
, mCosZenithRad( std::cos( static_cast<float>( lightAngle * M_PI ) / 180.0f ) )
27+
, mSinZenithRad( std::sin( static_cast<float>( lightAngle * M_PI ) / 180.0f ) )
28+
, mAzimuthRad( static_cast<float>( lightAzimuth * M_PI ) / 180.0f )
2629
{
2730
}
2831

2932
float QgsHillshadeFilter::processNineCellWindow( float *x11, float *x21, float *x31,
3033
float *x12, float *x22, float *x32,
3134
float *x13, float *x23, float *x33 )
3235
{
36+
3337
float derX = calcFirstDerX( x11, x21, x31, x12, x22, x32, x13, x23, x33 );
3438
float derY = calcFirstDerY( x11, x21, x31, x12, x22, x32, x13, x23, x33 );
3539

@@ -38,17 +42,43 @@ float QgsHillshadeFilter::processNineCellWindow( float *x11, float *x21, float *
3842
return mOutputNodataValue;
3943
}
4044

41-
float zenith_rad = mLightAngle * M_PI / 180.0;
4245
float slope_rad = std::atan( std::sqrt( derX * derX + derY * derY ) );
43-
float azimuth_rad = mLightAzimuth * M_PI / 180.0;
4446
float aspect_rad = 0;
4547
if ( derX == 0 && derY == 0 ) //aspect undefined, take a neutral value. Better solutions?
4648
{
47-
aspect_rad = azimuth_rad / 2.0;
49+
aspect_rad = mAzimuthRad / 2.0f;
4850
}
4951
else
5052
{
5153
aspect_rad = M_PI + std::atan2( derX, derY );
5254
}
53-
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 ) ) ) );
55+
return std::max( 0.0f, 255.0f * ( ( mCosZenithRad * std::cos( slope_rad ) ) +
56+
( mSinZenithRad * std::sin( slope_rad ) *
57+
std::cos( mAzimuthRad - aspect_rad ) ) ) );
58+
}
59+
60+
void QgsHillshadeFilter::setLightAzimuth( float azimuth )
61+
{
62+
mLightAzimuth = azimuth;
63+
mAzimuthRad = azimuth * static_cast<float>( M_PI ) / 180.0f;
64+
}
65+
66+
void QgsHillshadeFilter::setLightAngle( float angle )
67+
{
68+
mLightAngle = angle;
69+
mCosZenithRad = std::cos( angle * static_cast<float>( M_PI ) / 180.0f );
70+
mSinZenithRad = std::sin( angle * static_cast<float>( M_PI ) / 180.0f );
5471
}
72+
73+
#ifdef HAVE_OPENCL
74+
75+
void QgsHillshadeFilter::addExtraRasterParams( std::vector<float> &params )
76+
{
77+
78+
params.push_back( mCosZenithRad ); // cos_zenith_rad 5
79+
params.push_back( mSinZenithRad ); // sin_zenith_rad 6
80+
params.push_back( mAzimuthRad ); // azimuth_rad 7
81+
82+
}
83+
84+
#endif

‎src/analysis/raster/qgshillshadefilter.h

Lines changed: 22 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,33 @@ class ANALYSIS_EXPORT QgsHillshadeFilter: public QgsDerivativeFilter
3939
float *x13, float *x23, float *x33 ) override;
4040

4141
float lightAzimuth() const { return mLightAzimuth; }
42-
void setLightAzimuth( float azimuth ) { mLightAzimuth = azimuth; }
42+
void setLightAzimuth( float azimuth );
4343
float lightAngle() const { return mLightAngle; }
44-
void setLightAngle( float angle ) { mLightAngle = angle; }
44+
void setLightAngle( float angle );
4545

4646
private:
47+
48+
#ifdef HAVE_OPENCL
49+
50+
const QString openClProgramBaseName() const override
51+
{
52+
return QStringLiteral( "hillshade" );
53+
}
54+
#endif
55+
4756
float mLightAzimuth;
4857
float mLightAngle;
58+
// Precalculate for speed:
59+
float mCosZenithRad;
60+
float mSinZenithRad;
61+
float mAzimuthRad;
62+
63+
64+
#ifdef HAVE_OPENCL
65+
private:
66+
67+
void addExtraRasterParams( std::vector<float> &params ) override;
68+
#endif
4969

5070
};
5171

‎src/analysis/raster/qgsninecellfilter.cpp

Lines changed: 21 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -21,12 +21,18 @@
2121
#include "qgsfeedback.h"
2222
#include "qgsogrutils.h"
2323
#include "qgsmessagelog.h"
24+
25+
#ifdef HAVE_OPENCL
26+
#include "qgsopenclutils.h"
27+
#endif
28+
2429
#include <QFile>
2530
#include <QDebug>
2631
#include <QFileInfo>
2732
#include <iterator>
2833

2934

35+
3036
QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
3137
: mInputFile( inputFile )
3238
, mOutputFile( outputFile )
@@ -157,6 +163,7 @@ gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDa
157163
return outputDataset;
158164
}
159165

166+
#ifdef HAVE_OPENCL
160167

161168
int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *feedback )
162169
{
@@ -296,7 +303,8 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
296303
}
297304
queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
298305
}
299-
else // Overwrite from input, skip first and last
306+
else // Read line i + 1 and put it into scanline 3
307+
// Overwrite from input, skip first and last
300308
{
301309
if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
302310
{
@@ -334,7 +342,7 @@ int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *fee
334342
}
335343
return 0;
336344
}
337-
345+
#endif
338346

339347
int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
340348
{
@@ -393,7 +401,7 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
393401
float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
394402

395403
//values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
396-
for ( int i = 0; i < ySize; ++i )
404+
for ( int yIndex = 0; yIndex < ySize; ++yIndex )
397405
{
398406
if ( feedback && feedback->isCanceled() )
399407
{
@@ -402,10 +410,10 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
402410

403411
if ( feedback )
404412
{
405-
feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize );
413+
feedback->setProgress( 100.0 * static_cast< double >( yIndex ) / ySize );
406414
}
407415

408-
if ( i == 0 )
416+
if ( yIndex == 0 )
409417
{
410418
//fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
411419
for ( int a = 0; a < xSize + 2 ; ++a )
@@ -428,7 +436,7 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
428436
}
429437

430438
// Read scanline 3
431-
if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
439+
if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values
432440
{
433441
for ( int a = 0; a < xSize + 2; ++a )
434442
{
@@ -437,7 +445,7 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
437445
}
438446
else
439447
{
440-
if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
448+
if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
441449
{
442450
QgsDebugMsg( "Raster IO Error" );
443451
}
@@ -450,15 +458,16 @@ int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
450458

451459

452460
// j is the x axis index, skip 0 and last cell that have been filled with nodata
453-
for ( int j = 0; j < xSize ; ++j )
461+
for ( int xIndex = 0; xIndex < xSize ; ++xIndex )
454462
{
455-
resultLine[ j ] = processNineCellWindow( &scanLine1[ j ], &scanLine1[ j + 1 ], &scanLine1[ j + 2 ],
456-
&scanLine2[ j ], &scanLine2[ j + 1 ], &scanLine2[ j + 2 ],
457-
&scanLine3[ j ], &scanLine3[ j + 1 ], &scanLine3[ j + 2 ] );
463+
// cells(x, y) x11, x21, x31, x12, x22, x32, x13, x23, x33
464+
resultLine[ xIndex ] = processNineCellWindow( &scanLine1[ xIndex ], &scanLine1[ xIndex + 1 ], &scanLine1[ xIndex + 2 ],
465+
&scanLine2[ xIndex ], &scanLine2[ xIndex + 1 ], &scanLine2[ xIndex + 2 ],
466+
&scanLine3[ xIndex ], &scanLine3[ xIndex + 1 ], &scanLine3[ xIndex + 2 ] );
458467

459468
}
460469

461-
if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
470+
if ( GDALRasterIO( outputRasterBand, GF_Write, 0, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
462471
{
463472
QgsDebugMsg( "Raster IO Error" );
464473
}

‎src/analysis/raster/qgsninecellfilter.h

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,6 @@
2222
#include "gdal.h"
2323
#include "qgis_analysis.h"
2424
#include "qgsogrutils.h"
25-
#include "qgsopenclutils.h"
2625

2726
class QgsFeedback;
2827

‎src/analysis/raster/qgsruggednessfilter.cpp

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -24,11 +24,6 @@ QgsRuggednessFilter::QgsRuggednessFilter( const QString &inputFile, const QStrin
2424

2525
}
2626

27-
QgsRuggednessFilter::QgsRuggednessFilter()
28-
: QgsNineCellFilter( QString(), QString(), QString() )
29-
{
30-
31-
}
3227

3328
float QgsRuggednessFilter::processNineCellWindow( float *x11, float *x21, float *x31,
3429
float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 )

‎src/analysis/raster/qgsruggednessfilter.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,15 @@ class ANALYSIS_EXPORT QgsRuggednessFilter: public QgsNineCellFilter
3939
float *x12, float *x22, float *x32,
4040
float *x13, float *x23, float *x33 ) override;
4141

42+
#ifdef HAVE_OPENCL
4243
private:
4344
QgsRuggednessFilter();
4445

4546
virtual const QString openClProgramBaseName() const override
4647
{
4748
return QStringLiteral( "ruggedness" );
4849
}
50+
#endif
4951

5052
};
5153

‎src/analysis/raster/qgsslopefilter.h

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,12 +36,15 @@ class ANALYSIS_EXPORT QgsSlopeFilter: public QgsDerivativeFilter
3636
float *x12, float *x22, float *x32,
3737
float *x13, float *x23, float *x33 ) override;
3838

39+
40+
#ifdef HAVE_OPENCL
3941
private:
4042

4143
virtual const QString openClProgramBaseName() const override
4244
{
4345
return QStringLiteral( "slope" );
4446
}
47+
#endif
4548

4649
};
4750

0 commit comments

Comments
 (0)
Please sign in to comment.