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