QGIS API Documentation 3.41.0-Master (d5b93354e9c)
Loading...
Searching...
No Matches
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
37QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
38 : mInputFile( inputFile )
39 , mOutputFile( outputFile )
40 , mOutputFormat( outputFormat )
41{
42}
43
44// TODO: return an anum instead of an int
46{
47#ifdef HAVE_OPENCL
48 if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && !openClProgramBaseName().isEmpty() )
49 {
50 // Load the program sources
51 const QString source( QgsOpenClUtils::sourceFromBaseName( openClProgramBaseName() ) );
52 if ( !source.isEmpty() )
53 {
54 try
55 {
56 QgsDebugMsgLevel( QStringLiteral( "Running OpenCL program: %1" ).arg( openClProgramBaseName() ), 2 );
57 return processRasterGPU( source, feedback );
58 }
59 catch ( cl::Error &e )
60 {
61 const QString err = QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what(), QgsOpenClUtils::errorText( e.err() ) );
63 throw QgsProcessingException( err );
64 }
65 }
66 else
67 {
68 const QString err = QObject::tr( "Error loading OpenCL program sources" );
70 throw QgsProcessingException( err );
71 }
72 }
73 else
74 {
75 return processRasterCPU( feedback );
76 }
77#ifndef _MSC_VER
78 return 1;
79#endif
80#else
81 return processRasterCPU( feedback );
82#endif
83}
84
85gdal::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
102GDALDriverH 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
121gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
122{
123 if ( !inputDataset )
124 {
125 return nullptr;
126 }
127
128 const int xSize = GDALGetRasterXSize( inputDataset );
129 const 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 {
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
168int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *feedback )
169{
170 GDALAllRegister();
171
172 //open input file
173 int xSize, ySize;
174 const gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
175 if ( !inputDataset )
176 {
177 return 1; //opening of input file failed
178 }
179
180 //output driver
181 GDALDriverH outputDriver = openOutputDriver();
182 if ( !outputDriver )
183 {
184 return 2;
185 }
186
187 gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
188 if ( !outputDataset )
189 {
190 return 3; //create operation on output file failed
191 }
192
193 //open first raster band for reading (operation is only for single band raster)
194 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
195 if ( !rasterBand )
196 {
197 return 4;
198 }
199 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
200
201 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
202 if ( !outputRasterBand )
203 {
204 return 5;
205 }
206 //try to set -9999 as nodata value
207 GDALSetRasterNoDataValue( outputRasterBand, -9999 );
208 mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
209
210 if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
211 {
212 return 6;
213 }
214
215 // Prepare context and queue
216 const cl::Context ctx = QgsOpenClUtils::context();
217 cl::CommandQueue queue = QgsOpenClUtils::commandQueue();
218
219 //keep only three scanlines in memory at a time, make room for initial and final nodata
220 QgsOpenClUtils::CPLAllocator<float> scanLine( xSize + 2 );
221 QgsOpenClUtils::CPLAllocator<float> resultLine( xSize );
222
223 // Cast to float (because double just crashes on some GPUs)
224 std::vector<float> rasterParams;
225
226 rasterParams.push_back( mInputNodataValue ); // 0
227 rasterParams.push_back( mOutputNodataValue ); // 1
228 rasterParams.push_back( mZFactor ); // 2
229 rasterParams.push_back( mCellSizeX ); // 3
230 rasterParams.push_back( mCellSizeY ); // 4
231
232 // Allow subclasses to add extra params needed for computation:
233 // used to pass additional args to opencl program
234 addExtraRasterParams( rasterParams );
235
236 const std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
237 const std::size_t inputSize( sizeof( float ) * ( xSize ) );
238
239 cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr );
240 cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
241 cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
242 cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
243 cl::Buffer *scanLineBuffer[3] = { &scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer };
244 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, inputSize, nullptr, nullptr );
245
246 // Create a program from the kernel source
247 const cl::Program program( QgsOpenClUtils::buildProgram( source, QgsOpenClUtils::ExceptionBehavior::Throw ) );
248
249 // Create the OpenCL kernel
250 auto kernel = cl::KernelFunctor<
251 cl::Buffer &,
252 cl::Buffer &,
253 cl::Buffer &,
254 cl::Buffer &,
255 cl::Buffer &>( program, "processNineCellWindow" );
256
257 // Rotate buffer index
258 std::vector<int> rowIndex = { 0, 1, 2 };
259
260 // values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
261 for ( int i = 0; i < ySize; ++i )
262 {
263 if ( feedback && feedback->isCanceled() )
264 {
265 break;
266 }
267
268 if ( feedback )
269 {
270 feedback->setProgress( 100.0 * static_cast<double>( i ) / ySize );
271 }
272
273 if ( i == 0 )
274 {
275 // Fill scanline 1 with (input) nodata for the values above the first row and
276 // feed scanline2 with the first actual data row
277 for ( int a = 0; a < xSize + 2; ++a )
278 {
279 scanLine[a] = mInputNodataValue;
280 }
281 queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
282
283 // Read scanline2: first real raster row
284 if ( GDALRasterIO( rasterBand, GF_Read, 0, i, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
285 {
286 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
287 }
288 queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
289
290 // Read scanline3: second real raster row
291 if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
292 {
293 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
294 }
295 queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
296 }
297 else
298 {
299 // Normally fetch only scanLine3 and move forward one row
300 // Read scanline 3, fill the last row with nodata values if it's the last iteration
301 if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
302 {
303 for ( int a = 0; a < xSize + 2; ++a )
304 {
305 scanLine[a] = mInputNodataValue;
306 }
307 queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
308 }
309 else // Read line i + 1 and put it into scanline 3
310 // Overwrite from input, skip first and last
311 {
312 if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
313 {
314 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
315 }
316 queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
317 }
318 }
319
320 kernel( cl::EnqueueArgs( queue, cl::NDRange( xSize ) ), *scanLineBuffer[rowIndex[0]], *scanLineBuffer[rowIndex[1]], *scanLineBuffer[rowIndex[2]], resultLineBuffer, rasterParamsBuffer );
321
322 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, inputSize, resultLine.get() );
323
324 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine.get(), xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
325 {
326 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
327 }
328 std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
329 }
330
331 if ( feedback && feedback->isCanceled() )
332 {
333 //delete the dataset without closing (because it is faster)
334 gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
335 return 7;
336 }
337 return 0;
338}
339#endif
340
341
342// TODO: return an anum instead of an int
343int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
344{
345 GDALAllRegister();
346
347 //open input file
348 int xSize, ySize;
349 const gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
350 if ( !inputDataset )
351 {
352 return 1; //opening of input file failed
353 }
354
355 //output driver
356 GDALDriverH outputDriver = openOutputDriver();
357 if ( !outputDriver )
358 {
359 return 2;
360 }
361
362 gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
363 if ( !outputDataset )
364 {
365 return 3; //create operation on output file failed
366 }
367
368 //open first raster band for reading (operation is only for single band raster)
369 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
370 if ( !rasterBand )
371 {
372 return 4;
373 }
374 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
375
376 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
377 if ( !outputRasterBand )
378 {
379 return 5;
380 }
381 //try to set -9999 as nodata value
382 GDALSetRasterNoDataValue( outputRasterBand, -9999 );
383 mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
384
385 if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
386 {
387 return 6;
388 }
389
390 //keep only three scanlines in memory at a time, make room for initial and final nodata
391 const std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
392 float *scanLine1 = ( float * ) CPLMalloc( bufferSize );
393 float *scanLine2 = ( float * ) CPLMalloc( bufferSize );
394 float *scanLine3 = ( float * ) CPLMalloc( bufferSize );
395
396 float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
397
398 //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
399 for ( int yIndex = 0; yIndex < ySize; ++yIndex )
400 {
401 if ( feedback && feedback->isCanceled() )
402 {
403 break;
404 }
405
406 if ( feedback )
407 {
408 feedback->setProgress( 100.0 * static_cast<double>( yIndex ) / ySize );
409 }
410
411 if ( yIndex == 0 )
412 {
413 //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
414 for ( int a = 0; a < xSize + 2; ++a )
415 {
416 scanLine1[a] = mInputNodataValue;
417 }
418 // Read scanline2
419 if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, &scanLine2[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
420 {
421 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
422 }
423 }
424 else
425 {
426 //normally fetch only scanLine3 and release scanline 1 if we move forward one row
427 CPLFree( scanLine1 );
428 scanLine1 = scanLine2;
429 scanLine2 = scanLine3;
430 scanLine3 = ( float * ) CPLMalloc( bufferSize );
431 }
432
433 // Read scanline 3
434 if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values
435 {
436 for ( int a = 0; a < xSize + 2; ++a )
437 {
438 scanLine3[a] = mInputNodataValue;
439 }
440 }
441 else
442 {
443 if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
444 {
445 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
446 }
447 }
448 // Set first and last extra columns to nodata
449 scanLine1[0] = scanLine1[xSize + 1] = mInputNodataValue;
450 scanLine2[0] = scanLine2[xSize + 1] = mInputNodataValue;
451 scanLine3[0] = scanLine3[xSize + 1] = mInputNodataValue;
452
453
454 // j is the x axis index, skip 0 and last cell that have been filled with nodata
455 for ( int xIndex = 0; xIndex < xSize; ++xIndex )
456 {
457 // cells(x, y) x11, x21, x31, x12, x22, x32, x13, x23, x33
458 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] );
459 }
460
461 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
462 {
463 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
464 }
465 }
466
467 CPLFree( resultLine );
468 CPLFree( scanLine1 );
469 CPLFree( scanLine2 );
470 CPLFree( scanLine3 );
471
472 if ( feedback && feedback->isCanceled() )
473 {
474 //delete the dataset without closing (because it is faster)
475 gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
476 return 7;
477 }
478 return 0;
479}
@ Critical
Critical/error message.
Definition qgis.h:157
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
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 ...
@ Throw
Write errors in the message log and re-throw exceptions.
static QLatin1String LOGMESSAGE_TAG
OpenCL string for message logs.
Custom exception class for processing related exceptions.
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.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
void * GDALDatasetH
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39
#define QgsDebugError(str)
Definition qgslogger.h:38
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...