Branch data Line data Source code
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 : 0 : QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
38 : 0 : : mInputFile( inputFile )
39 : 0 : , mOutputFile( outputFile )
40 : 0 : , mOutputFormat( outputFormat )
41 : 0 : {
42 : :
43 : 0 : }
44 : :
45 : : // TODO: return an anum instead of an int
46 : 0 : 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 : 0 : return processRasterCPU( feedback );
84 : : #endif
85 : : }
86 : :
87 : 0 : gdal::dataset_unique_ptr QgsNineCellFilter::openInputFile( int &nCellsX, int &nCellsY )
88 : : {
89 : 0 : gdal::dataset_unique_ptr inputDataset( GDALOpen( mInputFile.toUtf8().constData(), GA_ReadOnly ) );
90 : 0 : if ( inputDataset )
91 : : {
92 : 0 : nCellsX = GDALGetRasterXSize( inputDataset.get() );
93 : 0 : nCellsY = GDALGetRasterYSize( inputDataset.get() );
94 : :
95 : : //we need at least one band
96 : 0 : if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
97 : : {
98 : 0 : return nullptr;
99 : : }
100 : 0 : }
101 : 0 : return inputDataset;
102 : 0 : }
103 : :
104 : 0 : GDALDriverH QgsNineCellFilter::openOutputDriver()
105 : : {
106 : : //open driver
107 : 0 : GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
108 : :
109 : 0 : if ( !outputDriver )
110 : : {
111 : 0 : return outputDriver; //return nullptr, driver does not exist
112 : : }
113 : :
114 : 0 : if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
115 : : {
116 : 0 : return nullptr; //driver exist, but it does not support the create operation
117 : : }
118 : :
119 : 0 : return outputDriver;
120 : 0 : }
121 : :
122 : :
123 : 0 : gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
124 : : {
125 : 0 : if ( !inputDataset )
126 : : {
127 : 0 : return nullptr;
128 : : }
129 : :
130 : 0 : int xSize = GDALGetRasterXSize( inputDataset );
131 : 0 : int ySize = GDALGetRasterYSize( inputDataset );
132 : :
133 : : //open output file
134 : 0 : char **papszOptions = nullptr;
135 : 0 : gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 1, GDT_Float32, papszOptions ) );
136 : 0 : if ( !outputDataset )
137 : : {
138 : 0 : return outputDataset;
139 : : }
140 : :
141 : : //get geotransform from inputDataset
142 : : double geotransform[6];
143 : 0 : if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
144 : 0 : {
145 : 0 : return nullptr;
146 : 0 : }
147 : 0 : GDALSetGeoTransform( outputDataset.get(), geotransform );
148 : 0 :
149 : : //make sure mCellSizeX and mCellSizeY are always > 0
150 : 0 : mCellSizeX = geotransform[1];
151 : 0 : if ( mCellSizeX < 0 )
152 : : {
153 : 0 : mCellSizeX = -mCellSizeX;
154 : 0 : }
155 : 0 : mCellSizeY = geotransform[5];
156 : 0 : if ( mCellSizeY < 0 )
157 : : {
158 : 0 : mCellSizeY = -mCellSizeY;
159 : 0 : }
160 : :
161 : 0 : const char *projection = GDALGetProjectionRef( inputDataset );
162 : 0 : GDALSetProjection( outputDataset.get(), projection );
163 : :
164 : 0 : return outputDataset;
165 : 0 : }
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 : 0 : int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
357 : : {
358 : :
359 : 0 : GDALAllRegister();
360 : :
361 : : //open input file
362 : : int xSize, ySize;
363 : 0 : gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
364 : 0 : if ( !inputDataset )
365 : : {
366 : 0 : return 1; //opening of input file failed
367 : : }
368 : :
369 : : //output driver
370 : 0 : GDALDriverH outputDriver = openOutputDriver();
371 : 0 : if ( !outputDriver )
372 : : {
373 : 0 : return 2;
374 : : }
375 : :
376 : 0 : gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
377 : 0 : if ( !outputDataset )
378 : : {
379 : 0 : 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 : 0 : GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
384 : 0 : if ( !rasterBand )
385 : : {
386 : 0 : return 4;
387 : : }
388 : 0 : mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
389 : :
390 : 0 : GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
391 : 0 : if ( !outputRasterBand )
392 : : {
393 : 0 : return 5;
394 : : }
395 : : //try to set -9999 as nodata value
396 : 0 : GDALSetRasterNoDataValue( outputRasterBand, -9999 );
397 : 0 : mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
398 : :
399 : 0 : if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
400 : : {
401 : 0 : return 6;
402 : : }
403 : :
404 : : //keep only three scanlines in memory at a time, make room for initial and final nodata
405 : 0 : std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
406 : 0 : float *scanLine1 = ( float * ) CPLMalloc( bufferSize );
407 : 0 : float *scanLine2 = ( float * ) CPLMalloc( bufferSize );
408 : 0 : float *scanLine3 = ( float * ) CPLMalloc( bufferSize );
409 : :
410 : 0 : 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 : 0 : for ( int yIndex = 0; yIndex < ySize; ++yIndex )
414 : : {
415 : 0 : if ( feedback && feedback->isCanceled() )
416 : : {
417 : 0 : break;
418 : : }
419 : :
420 : 0 : if ( feedback )
421 : : {
422 : 0 : feedback->setProgress( 100.0 * static_cast< double >( yIndex ) / ySize );
423 : 0 : }
424 : :
425 : 0 : 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 : 0 : for ( int a = 0; a < xSize + 2 ; ++a )
429 : : {
430 : 0 : scanLine1[a] = mInputNodataValue;
431 : 0 : }
432 : : // Read scanline2
433 : 0 : if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, &scanLine2[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
434 : : {
435 : 0 : QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
436 : 0 : }
437 : 0 : }
438 : : else
439 : : {
440 : : //normally fetch only scanLine3 and release scanline 1 if we move forward one row
441 : 0 : CPLFree( scanLine1 );
442 : 0 : scanLine1 = scanLine2;
443 : 0 : scanLine2 = scanLine3;
444 : 0 : scanLine3 = ( float * ) CPLMalloc( bufferSize );
445 : : }
446 : :
447 : : // Read scanline 3
448 : 0 : if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values
449 : : {
450 : 0 : for ( int a = 0; a < xSize + 2; ++a )
451 : : {
452 : 0 : scanLine3[a] = mInputNodataValue;
453 : 0 : }
454 : 0 : }
455 : : else
456 : : {
457 : 0 : if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
458 : : {
459 : 0 : QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
460 : 0 : }
461 : : }
462 : : // Set first and last extra columns to nodata
463 : 0 : scanLine1[0] = scanLine1[xSize + 1] = mInputNodataValue;
464 : 0 : scanLine2[0] = scanLine2[xSize + 1] = mInputNodataValue;
465 : 0 : 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 : 0 : for ( int xIndex = 0; xIndex < xSize ; ++xIndex )
471 : : {
472 : : // cells(x, y) x11, x21, x31, x12, x22, x32, x13, x23, x33
473 : 0 : resultLine[ xIndex ] = processNineCellWindow( &scanLine1[ xIndex ], &scanLine1[ xIndex + 1 ], &scanLine1[ xIndex + 2 ],
474 : 0 : &scanLine2[ xIndex ], &scanLine2[ xIndex + 1 ], &scanLine2[ xIndex + 2 ],
475 : 0 : &scanLine3[ xIndex ], &scanLine3[ xIndex + 1 ], &scanLine3[ xIndex + 2 ] );
476 : :
477 : 0 : }
478 : :
479 : 0 : if ( GDALRasterIO( outputRasterBand, GF_Write, 0, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
480 : : {
481 : 0 : QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
482 : 0 : }
483 : 0 : }
484 : :
485 : 0 : CPLFree( resultLine );
486 : 0 : CPLFree( scanLine1 );
487 : 0 : CPLFree( scanLine2 );
488 : 0 : CPLFree( scanLine3 );
489 : :
490 : 0 : if ( feedback && feedback->isCanceled() )
491 : : {
492 : : //delete the dataset without closing (because it is faster)
493 : 0 : gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
494 : 0 : return 7;
495 : : }
496 : 0 : return 0;
497 : 0 : }
|
