QGIS API Documentation  3.14.0-Pi (9f7028fd23)
qgsninecellfilter.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsninecellfilter.h - description
3  -------------------
4  begin : August 6th, 2009
5  copyright : (C) 2009 by Marco Hugentobler
6  email : marco dot hugentobler at karto dot baug dot ethz dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgsgdalutils.h"
19 #include "qgsninecellfilter.h"
20 #include "qgslogger.h"
21 #include "cpl_string.h"
22 #include "qgsfeedback.h"
23 #include "qgsogrutils.h"
24 #include "qgsmessagelog.h"
25 
26 #ifdef HAVE_OPENCL
27 #include "qgsopenclutils.h"
28 #endif
29 
30 #include <QFile>
31 #include <QDebug>
32 #include <QFileInfo>
33 #include <iterator>
34 
35 
36 
37 QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
38  : mInputFile( inputFile )
39  , mOutputFile( outputFile )
40  , mOutputFormat( outputFormat )
41 {
42 
43 }
44 
45 // TODO: return an anum instead of an int
47 {
48 #ifdef HAVE_OPENCL
49  if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! openClProgramBaseName( ).isEmpty() )
50  {
51  // Load the program sources
52  QString source( QgsOpenClUtils::sourceFromBaseName( openClProgramBaseName( ) ) );
53  if ( ! source.isEmpty() )
54  {
55  try
56  {
57  QgsDebugMsg( QObject::tr( "Running OpenCL program: %1" ).arg( openClProgramBaseName( ) ) );
58  return processRasterGPU( source, feedback );
59  }
60  catch ( cl::Error &e )
61  {
62  QString err = QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what( ), QgsOpenClUtils::errorText( e.err( ) ) );
64  throw QgsProcessingException( err );
65  }
66  }
67  else
68  {
69  QString err = QObject::tr( "Error loading OpenCL program sources" );
72  throw QgsProcessingException( err );
73  }
74  }
75  else
76  {
77  return processRasterCPU( feedback );
78  }
79 #ifndef _MSC_VER
80  return 1;
81 #endif
82 #else
83  return processRasterCPU( feedback );
84 #endif
85 }
86 
87 gdal::dataset_unique_ptr QgsNineCellFilter::openInputFile( int &nCellsX, int &nCellsY )
88 {
89  gdal::dataset_unique_ptr inputDataset( GDALOpen( mInputFile.toUtf8().constData(), GA_ReadOnly ) );
90  if ( inputDataset )
91  {
92  nCellsX = GDALGetRasterXSize( inputDataset.get() );
93  nCellsY = GDALGetRasterYSize( inputDataset.get() );
94 
95  //we need at least one band
96  if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
97  {
98  return nullptr;
99  }
100  }
101  return inputDataset;
102 }
103 
104 GDALDriverH QgsNineCellFilter::openOutputDriver()
105 {
106  //open driver
107  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
108 
109  if ( !outputDriver )
110  {
111  return outputDriver; //return nullptr, driver does not exist
112  }
113 
114  if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
115  {
116  return nullptr; //driver exist, but it does not support the create operation
117  }
118 
119  return outputDriver;
120 }
121 
122 
123 gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
124 {
125  if ( !inputDataset )
126  {
127  return nullptr;
128  }
129 
130  int xSize = GDALGetRasterXSize( inputDataset );
131  int ySize = GDALGetRasterYSize( inputDataset );
132 
133  //open output file
134  char **papszOptions = nullptr;
135  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 1, GDT_Float32, papszOptions ) );
136  if ( !outputDataset )
137  {
138  return outputDataset;
139  }
140 
141  //get geotransform from inputDataset
142  double geotransform[6];
143  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
144  {
145  return nullptr;
146  }
147  GDALSetGeoTransform( outputDataset.get(), geotransform );
148 
149  //make sure mCellSizeX and mCellSizeY are always > 0
150  mCellSizeX = geotransform[1];
151  if ( mCellSizeX < 0 )
152  {
154  }
155  mCellSizeY = geotransform[5];
156  if ( mCellSizeY < 0 )
157  {
159  }
160 
161  const char *projection = GDALGetProjectionRef( inputDataset );
162  GDALSetProjection( outputDataset.get(), projection );
163 
164  return outputDataset;
165 }
166 
167 #ifdef HAVE_OPENCL
168 
169 // TODO: return an anum instead of an int
170 int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *feedback )
171 {
172 
173  GDALAllRegister();
174 
175  //open input file
176  int xSize, ySize;
177  gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
178  if ( !inputDataset )
179  {
180  return 1; //opening of input file failed
181  }
182 
183  //output driver
184  GDALDriverH outputDriver = openOutputDriver();
185  if ( !outputDriver )
186  {
187  return 2;
188  }
189 
190  gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
191  if ( !outputDataset )
192  {
193  return 3; //create operation on output file failed
194  }
195 
196  //open first raster band for reading (operation is only for single band raster)
197  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
198  if ( !rasterBand )
199  {
200  return 4;
201  }
202  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
203 
204  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
205  if ( !outputRasterBand )
206  {
207  return 5;
208  }
209  //try to set -9999 as nodata value
210  GDALSetRasterNoDataValue( outputRasterBand, -9999 );
211  mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
212 
213  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
214  {
215  return 6;
216  }
217 
218  // Prepare context and queue
219  cl::Context ctx = QgsOpenClUtils::context();
220  cl::CommandQueue queue = QgsOpenClUtils::commandQueue();
221 
222  //keep only three scanlines in memory at a time, make room for initial and final nodata
223  QgsOpenClUtils::CPLAllocator<float> scanLine( xSize + 2 );
224  QgsOpenClUtils::CPLAllocator<float> resultLine( xSize );
225 
226  // Cast to float (because double just crashes on some GPUs)
227  std::vector<float> rasterParams;
228 
229  rasterParams.push_back( mInputNodataValue ); // 0
230  rasterParams.push_back( mOutputNodataValue ); // 1
231  rasterParams.push_back( mZFactor ); // 2
232  rasterParams.push_back( mCellSizeX ); // 3
233  rasterParams.push_back( mCellSizeY ); // 4
234 
235  // Allow subclasses to add extra params needed for computation:
236  // used to pass additional args to opencl program
237  addExtraRasterParams( rasterParams );
238 
239  std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
240  std::size_t inputSize( sizeof( float ) * ( xSize ) );
241 
242  cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr );
243  cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
244  cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
245  cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
246  cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer};
247  cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, inputSize, nullptr, nullptr );
248 
249  // Create a program from the kernel source
250  cl::Program program( QgsOpenClUtils::buildProgram( source, QgsOpenClUtils::ExceptionBehavior::Throw ) );
251 
252  // Create the OpenCL kernel
253  auto kernel = cl::KernelFunctor <
254  cl::Buffer &,
255  cl::Buffer &,
256  cl::Buffer &,
257  cl::Buffer &,
258  cl::Buffer &
259  > ( program, "processNineCellWindow" );
260 
261  // Rotate buffer index
262  std::vector<int> rowIndex = {0, 1, 2};
263 
264  // values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
265  for ( int i = 0; i < ySize; ++i )
266  {
267  if ( feedback && feedback->isCanceled() )
268  {
269  break;
270  }
271 
272  if ( feedback )
273  {
274  feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize );
275  }
276 
277  if ( i == 0 )
278  {
279  // Fill scanline 1 with (input) nodata for the values above the first row and
280  // feed scanline2 with the first actual data row
281  for ( int a = 0; a < xSize + 2 ; ++a )
282  {
283  scanLine[a] = mInputNodataValue;
284  }
285  queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
286 
287  // Read scanline2: first real raster row
288  if ( GDALRasterIO( rasterBand, GF_Read, 0, i, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
289  {
290  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
291  }
292  queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
293 
294  // Read scanline3: second real raster row
295  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
296  {
297  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
298  }
299  queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
300  }
301  else
302  {
303  // Normally fetch only scanLine3 and move forward one row
304  // Read scanline 3, fill the last row with nodata values if it's the last iteration
305  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
306  {
307  for ( int a = 0; a < xSize + 2; ++a )
308  {
309  scanLine[a] = mInputNodataValue;
310  }
311  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
312  }
313  else // Read line i + 1 and put it into scanline 3
314  // Overwrite from input, skip first and last
315  {
316  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
317  {
318  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
319  }
320  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
321  }
322  }
323 
324  kernel( cl::EnqueueArgs(
325  queue,
326  cl::NDRange( xSize )
327  ),
328  *scanLineBuffer[rowIndex[0]],
329  *scanLineBuffer[rowIndex[1]],
330  *scanLineBuffer[rowIndex[2]],
331  resultLineBuffer,
332  rasterParamsBuffer
333  );
334 
335  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, inputSize, resultLine.get() );
336 
337  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine.get(), xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
338  {
339  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
340  }
341  std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
342  }
343 
344  if ( feedback && feedback->isCanceled() )
345  {
346  //delete the dataset without closing (because it is faster)
347  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
348  return 7;
349  }
350  return 0;
351 }
352 #endif
353 
354 
355 // TODO: return an anum instead of an int
356 int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
357 {
358 
359  GDALAllRegister();
360 
361  //open input file
362  int xSize, ySize;
363  gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
364  if ( !inputDataset )
365  {
366  return 1; //opening of input file failed
367  }
368 
369  //output driver
370  GDALDriverH outputDriver = openOutputDriver();
371  if ( !outputDriver )
372  {
373  return 2;
374  }
375 
376  gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
377  if ( !outputDataset )
378  {
379  return 3; //create operation on output file failed
380  }
381 
382  //open first raster band for reading (operation is only for single band raster)
383  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
384  if ( !rasterBand )
385  {
386  return 4;
387  }
388  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
389 
390  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
391  if ( !outputRasterBand )
392  {
393  return 5;
394  }
395  //try to set -9999 as nodata value
396  GDALSetRasterNoDataValue( outputRasterBand, -9999 );
397  mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
398 
399  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
400  {
401  return 6;
402  }
403 
404  //keep only three scanlines in memory at a time, make room for initial and final nodata
405  std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
406  float *scanLine1 = ( float * ) CPLMalloc( bufferSize );
407  float *scanLine2 = ( float * ) CPLMalloc( bufferSize );
408  float *scanLine3 = ( float * ) CPLMalloc( bufferSize );
409 
410  float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
411 
412  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
413  for ( int yIndex = 0; yIndex < ySize; ++yIndex )
414  {
415  if ( feedback && feedback->isCanceled() )
416  {
417  break;
418  }
419 
420  if ( feedback )
421  {
422  feedback->setProgress( 100.0 * static_cast< double >( yIndex ) / ySize );
423  }
424 
425  if ( yIndex == 0 )
426  {
427  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
428  for ( int a = 0; a < xSize + 2 ; ++a )
429  {
430  scanLine1[a] = mInputNodataValue;
431  }
432  // Read scanline2
433  if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, &scanLine2[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
434  {
435  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
436  }
437  }
438  else
439  {
440  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
441  CPLFree( scanLine1 );
442  scanLine1 = scanLine2;
443  scanLine2 = scanLine3;
444  scanLine3 = ( float * ) CPLMalloc( bufferSize );
445  }
446 
447  // Read scanline 3
448  if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values
449  {
450  for ( int a = 0; a < xSize + 2; ++a )
451  {
452  scanLine3[a] = mInputNodataValue;
453  }
454  }
455  else
456  {
457  if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
458  {
459  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
460  }
461  }
462  // Set first and last extra columns to nodata
463  scanLine1[0] = scanLine1[xSize + 1] = mInputNodataValue;
464  scanLine2[0] = scanLine2[xSize + 1] = mInputNodataValue;
465  scanLine3[0] = scanLine3[xSize + 1] = mInputNodataValue;
466 
467 
468 
469  // j is the x axis index, skip 0 and last cell that have been filled with nodata
470  for ( int xIndex = 0; xIndex < xSize ; ++xIndex )
471  {
472  // cells(x, y) x11, x21, x31, x12, x22, x32, x13, x23, x33
473  resultLine[ xIndex ] = processNineCellWindow( &scanLine1[ xIndex ], &scanLine1[ xIndex + 1 ], &scanLine1[ xIndex + 2 ],
474  &scanLine2[ xIndex ], &scanLine2[ xIndex + 1 ], &scanLine2[ xIndex + 2 ],
475  &scanLine3[ xIndex ], &scanLine3[ xIndex + 1 ], &scanLine3[ xIndex + 2 ] );
476 
477  }
478 
479  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
480  {
481  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
482  }
483  }
484 
485  CPLFree( resultLine );
486  CPLFree( scanLine1 );
487  CPLFree( scanLine2 );
488  CPLFree( scanLine3 );
489 
490  if ( feedback && feedback->isCanceled() )
491  {
492  //delete the dataset without closing (because it is faster)
493  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
494  return 7;
495  }
496  return 0;
497 }
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:75
QgsOpenClUtils::errorText
static QString errorText(const int errorCode)
Returns a string representation from an OpenCL errorCode.
Definition: qgsopenclutils.cpp:510
qgsgdalutils.h
qgsogrutils.h
QgsGdalUtils::supportsRasterCreate
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
Definition: qgsgdalutils.cpp:30
QgsDebugMsg
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QgsNineCellFilter::processRaster
int processRaster(QgsFeedback *feedback=nullptr)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
Definition: qgsninecellfilter.cpp:46
QgsNineCellFilter::mInputFile
QString mInputFile
Definition: qgsninecellfilter.h:136
QgsOpenClUtils::available
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
Definition: qgsopenclutils.cpp:463
QgsNineCellFilter::mCellSizeX
double mCellSizeX
Definition: qgsninecellfilter.h:140
QgsOpenClUtils::CPLAllocator
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...
Definition: qgsopenclutils.h:237
gdal::dataset_unique_ptr
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:134
QgsFeedback
Definition: qgsfeedback.h:43
QgsNineCellFilter::processNineCellWindow
virtual float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33)=0
Calculates output value from nine input values.
QgsNineCellFilter::mOutputFormat
QString mOutputFormat
Definition: qgsninecellfilter.h:138
QgsOpenClUtils::sourceFromBaseName
static QString sourceFromBaseName(const QString &baseName)
Returns the full path to a an OpenCL source file from the baseName ('.cl' extension is automatically ...
Definition: qgsopenclutils.cpp:495
QgsOpenClUtils::enabled
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
Definition: qgsopenclutils.cpp:265
QgsNineCellFilter::mOutputFile
QString mOutputFile
Definition: qgsninecellfilter.h:137
QgsNineCellFilter::mOutputNodataValue
float mOutputNodataValue
The nodata value of the output layer.
Definition: qgsninecellfilter.h:145
gdal::fast_delete_and_close
void CORE_EXPORT fast_delete_and_close(dataset_unique_ptr &dataset, GDALDriverH driver, const QString &path)
Performs a fast close of an unwanted GDAL dataset handle by deleting the underlying data store.
Definition: qgsogrutils.cpp:71
QgsFeedback::isCanceled
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:66
QgsNineCellFilter::mInputNodataValue
float mInputNodataValue
The nodata value of the input layer.
Definition: qgsninecellfilter.h:143
QgsNineCellFilter::QgsNineCellFilter
QgsNineCellFilter(const QString &inputFile, const QString &outputFile, const QString &outputFormat)
Constructor that takes input file, output file and output format (GDAL string)
Definition: qgsninecellfilter.cpp:37
QgsNineCellFilter::mCellSizeY
double mCellSizeY
Definition: qgsninecellfilter.h:141
QgsNineCellFilter::mZFactor
double mZFactor
Scale factor for z-value if x-/y- units are different to z-units (111120 for degree->meters and 37040...
Definition: qgsninecellfilter.h:147
QgsMessageLog::logMessage
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
Definition: qgsmessagelog.cpp:27
QgsOpenClUtils::commandQueue
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
Definition: qgsopenclutils.cpp:607
QgsOpenClUtils::context
static cl::Context context()
Context factory.
Definition: qgsopenclutils.cpp:628
Qgis::Critical
@ Critical
Definition: qgis.h:105
qgsopenclutils.h
qgslogger.h
qgsfeedback.h
QgsOpenClUtils::LOGMESSAGE_TAG
static QLatin1String LOGMESSAGE_TAG
OpenCL string for message logs.
Definition: qgsopenclutils.h:183
QgsOpenClUtils::buildProgram
static Q_DECL_DEPRECATED cl::Program buildProgram(const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior=Catch)
Build the program from source in the given context and depending on exceptionBehavior can throw or ca...
Definition: qgsopenclutils.cpp:642
QgsProcessingException
Definition: qgsexception.h:82
qgsninecellfilter.h
qgsmessagelog.h
GDALDatasetH
void * GDALDatasetH
Definition: qgsalignraster.h:28