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