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
QgsNineCellFilter(const QString & inputFile,const QString & outputFile,const QString & outputFormat)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
processRaster(QgsFeedback * feedback)46 int QgsNineCellFilter::processRaster( QgsFeedback *feedback )
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( ) ) );
63 QgsMessageLog::logMessage( err, QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
64 throw QgsProcessingException( err );
65 }
66 }
67 else
68 {
69 QString err = QObject::tr( "Error loading OpenCL program sources" );
70 QgsMessageLog::logMessage( err,
71 QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::Critical );
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
openInputFile(int & nCellsX,int & nCellsY)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
openOutputDriver()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
openOutputFile(GDALDatasetH inputDataset,GDALDriverH outputDriver)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 {
153 mCellSizeX = -mCellSizeX;
154 }
155 mCellSizeY = geotransform[5];
156 if ( mCellSizeY < 0 )
157 {
158 mCellSizeY = -mCellSizeY;
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
processRasterGPU(const QString & source,QgsFeedback * feedback)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
processRasterCPU(QgsFeedback * feedback)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 }
498