Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsrasterface.cpp - Internal raster processing modules interface
3 : : --------------------------------------
4 : : Date : Jun 21, 2012
5 : : Copyright : (C) 2012 by Radim Blazek
6 : : email : radim dot blazek at gmail dot com
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 <limits>
19 : : #include <typeinfo>
20 : :
21 : : #include <QByteArray>
22 : : #include <QTime>
23 : : #include <QStringList>
24 : :
25 : : #include "qgslogger.h"
26 : : #include "qgsrasterbandstats.h"
27 : : #include "qgsrasterhistogram.h"
28 : : #include "qgsrasterinterface.h"
29 : : #include "qgsrectangle.h"
30 : :
31 : 0 : QgsRasterInterface::QgsRasterInterface( QgsRasterInterface *input )
32 : 0 : : mInput( input )
33 : 0 : {
34 : 0 : }
35 : :
36 : 0 : void QgsRasterInterface::initStatistics( QgsRasterBandStats &statistics,
37 : : int bandNo,
38 : : int stats,
39 : : const QgsRectangle &boundingBox,
40 : : int sampleSize )
41 : : {
42 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 sampleSize = %2" ).arg( bandNo ).arg( sampleSize ), 4 );
43 : :
44 : 0 : statistics.bandNumber = bandNo;
45 : 0 : statistics.statsGathered = stats;
46 : :
47 : 0 : QgsRectangle finalExtent;
48 : 0 : if ( boundingBox.isEmpty() )
49 : : {
50 : 0 : finalExtent = extent();
51 : 0 : }
52 : : else
53 : : {
54 : 0 : finalExtent = extent().intersect( boundingBox );
55 : : }
56 : 0 : statistics.extent = finalExtent;
57 : :
58 : 0 : if ( sampleSize > 0 )
59 : : {
60 : : // Calc resolution from theSampleSize
61 : : double xRes, yRes;
62 : 0 : xRes = yRes = std::sqrt( ( finalExtent.width() * finalExtent.height() ) / sampleSize );
63 : :
64 : : // But limit by physical resolution
65 : 0 : if ( capabilities() & Size )
66 : : {
67 : 0 : double srcXRes = extent().width() / xSize();
68 : 0 : double srcYRes = extent().height() / ySize();
69 : 0 : if ( xRes < srcXRes ) xRes = srcXRes;
70 : 0 : if ( yRes < srcYRes ) yRes = srcYRes;
71 : 0 : }
72 : 0 : QgsDebugMsgLevel( QStringLiteral( "xRes = %1 yRes = %2" ).arg( xRes ).arg( yRes ), 4 );
73 : :
74 : 0 : statistics.width = static_cast <int>( std::ceil( finalExtent.width() / xRes ) );
75 : 0 : statistics.height = static_cast <int>( std::ceil( finalExtent.height() / yRes ) );
76 : 0 : }
77 : : else
78 : : {
79 : 0 : if ( capabilities() & Size )
80 : : {
81 : 0 : statistics.width = xSize();
82 : 0 : statistics.height = ySize();
83 : 0 : }
84 : : else
85 : : {
86 : 0 : statistics.width = 1000;
87 : 0 : statistics.height = 1000;
88 : : }
89 : : }
90 : 0 : QgsDebugMsgLevel( QStringLiteral( "theStatistics.width = %1 statistics.height = %2" ).arg( statistics.width ).arg( statistics.height ), 4 );
91 : 0 : }
92 : :
93 : 0 : bool QgsRasterInterface::hasStatistics( int bandNo,
94 : : int stats,
95 : : const QgsRectangle &extent,
96 : : int sampleSize )
97 : : {
98 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 stats = %2 sampleSize = %3" ).arg( bandNo ).arg( stats ).arg( sampleSize ), 4 );
99 : 0 : if ( mStatistics.isEmpty() ) return false;
100 : :
101 : 0 : QgsRasterBandStats myRasterBandStats;
102 : 0 : initStatistics( myRasterBandStats, bandNo, stats, extent, sampleSize );
103 : :
104 : 0 : const auto constMStatistics = mStatistics;
105 : 0 : for ( const QgsRasterBandStats &stats : constMStatistics )
106 : : {
107 : 0 : if ( stats.contains( myRasterBandStats ) )
108 : : {
109 : 0 : QgsDebugMsgLevel( QStringLiteral( "Has cached statistics." ), 4 );
110 : 0 : return true;
111 : : }
112 : : }
113 : 0 : return false;
114 : 0 : }
115 : :
116 : 0 : QgsRasterBandStats QgsRasterInterface::bandStatistics( int bandNo,
117 : : int stats,
118 : : const QgsRectangle &extent,
119 : : int sampleSize, QgsRasterBlockFeedback *feedback )
120 : : {
121 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 stats = %2 sampleSize = %3" ).arg( bandNo ).arg( stats ).arg( sampleSize ), 4 );
122 : :
123 : : // TODO: null values set on raster layer!!!
124 : :
125 : 0 : QgsRasterBandStats myRasterBandStats;
126 : 0 : initStatistics( myRasterBandStats, bandNo, stats, extent, sampleSize );
127 : :
128 : 0 : const auto constMStatistics = mStatistics;
129 : 0 : for ( const QgsRasterBandStats &stats : constMStatistics )
130 : : {
131 : 0 : if ( stats.contains( myRasterBandStats ) )
132 : : {
133 : 0 : QgsDebugMsgLevel( QStringLiteral( "Using cached statistics." ), 4 );
134 : 0 : return stats;
135 : : }
136 : : }
137 : :
138 : 0 : QgsRectangle myExtent = myRasterBandStats.extent;
139 : 0 : int myWidth = myRasterBandStats.width;
140 : 0 : int myHeight = myRasterBandStats.height;
141 : :
142 : : //int myDataType = dataType( bandNo );
143 : :
144 : 0 : int myXBlockSize = xBlockSize();
145 : 0 : int myYBlockSize = yBlockSize();
146 : 0 : if ( myXBlockSize == 0 ) // should not happen, but happens
147 : : {
148 : 0 : myXBlockSize = 500;
149 : 0 : }
150 : 0 : if ( myYBlockSize == 0 ) // should not happen, but happens
151 : : {
152 : 0 : myYBlockSize = 500;
153 : 0 : }
154 : :
155 : 0 : int myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
156 : 0 : int myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
157 : :
158 : 0 : double myXRes = myExtent.width() / myWidth;
159 : 0 : double myYRes = myExtent.height() / myHeight;
160 : : // TODO: progress signals
161 : :
162 : : // used by single pass stdev
163 : 0 : double myMean = 0;
164 : 0 : double mySumOfSquares = 0;
165 : :
166 : 0 : bool myFirstIterationFlag = true;
167 : 0 : bool isNoData = false;
168 : 0 : for ( int myYBlock = 0; myYBlock < myNYBlocks; myYBlock++ )
169 : : {
170 : 0 : for ( int myXBlock = 0; myXBlock < myNXBlocks; myXBlock++ )
171 : : {
172 : 0 : if ( feedback && feedback->isCanceled() )
173 : 0 : return myRasterBandStats;
174 : :
175 : 0 : QgsDebugMsgLevel( QStringLiteral( "myYBlock = %1 myXBlock = %2" ).arg( myYBlock ).arg( myXBlock ), 4 );
176 : 0 : int myBlockWidth = std::min( myXBlockSize, myWidth - myXBlock * myXBlockSize );
177 : 0 : int myBlockHeight = std::min( myYBlockSize, myHeight - myYBlock * myYBlockSize );
178 : :
179 : 0 : double xmin = myExtent.xMinimum() + myXBlock * myXBlockSize * myXRes;
180 : 0 : double xmax = xmin + myBlockWidth * myXRes;
181 : 0 : double ymin = myExtent.yMaximum() - myYBlock * myYBlockSize * myYRes;
182 : 0 : double ymax = ymin - myBlockHeight * myYRes;
183 : :
184 : 0 : QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
185 : :
186 : 0 : std::unique_ptr< QgsRasterBlock > blk( block( bandNo, myPartExtent, myBlockWidth, myBlockHeight, feedback ) );
187 : :
188 : : // Collect the histogram counts.
189 : 0 : for ( qgssize i = 0; i < ( static_cast< qgssize >( myBlockHeight ) ) * myBlockWidth; i++ )
190 : : {
191 : 0 : double myValue = blk->valueAndNoData( i, isNoData );
192 : 0 : if ( isNoData )
193 : 0 : continue; // NULL
194 : :
195 : 0 : myRasterBandStats.sum += myValue;
196 : 0 : myRasterBandStats.elementCount++;
197 : :
198 : 0 : if ( !std::isfinite( myValue ) ) continue; // inf
199 : :
200 : 0 : if ( myFirstIterationFlag )
201 : : {
202 : 0 : myFirstIterationFlag = false;
203 : 0 : myRasterBandStats.minimumValue = myValue;
204 : 0 : myRasterBandStats.maximumValue = myValue;
205 : 0 : }
206 : : else
207 : : {
208 : 0 : if ( myValue < myRasterBandStats.minimumValue )
209 : : {
210 : 0 : myRasterBandStats.minimumValue = myValue;
211 : 0 : }
212 : 0 : if ( myValue > myRasterBandStats.maximumValue )
213 : : {
214 : 0 : myRasterBandStats.maximumValue = myValue;
215 : 0 : }
216 : : }
217 : :
218 : : // Single pass stdev
219 : 0 : double myDelta = myValue - myMean;
220 : 0 : myMean += myDelta / myRasterBandStats.elementCount;
221 : 0 : mySumOfSquares += myDelta * ( myValue - myMean );
222 : 0 : }
223 : 0 : }
224 : 0 : }
225 : :
226 : 0 : myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue;
227 : 0 : myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount;
228 : :
229 : 0 : myRasterBandStats.sumOfSquares = mySumOfSquares; // OK with single pass?
230 : :
231 : : // stdDev may differ from GDAL stats, because GDAL is using naive single pass
232 : : // algorithm which is more error prone (because of rounding errors)
233 : : // Divide result by sample size - 1 and get square root to get stdev
234 : 0 : myRasterBandStats.stdDev = std::sqrt( mySumOfSquares / ( myRasterBandStats.elementCount - 1 ) );
235 : :
236 : 0 : QgsDebugMsgLevel( QStringLiteral( "************ STATS **************" ), 4 );
237 : 0 : QgsDebugMsgLevel( QStringLiteral( "MIN %1" ).arg( myRasterBandStats.minimumValue ), 4 );
238 : 0 : QgsDebugMsgLevel( QStringLiteral( "MAX %1" ).arg( myRasterBandStats.maximumValue ), 4 );
239 : 0 : QgsDebugMsgLevel( QStringLiteral( "RANGE %1" ).arg( myRasterBandStats.range ), 4 );
240 : 0 : QgsDebugMsgLevel( QStringLiteral( "MEAN %1" ).arg( myRasterBandStats.mean ), 4 );
241 : 0 : QgsDebugMsgLevel( QStringLiteral( "STDDEV %1" ).arg( myRasterBandStats.stdDev ), 4 );
242 : :
243 : 0 : myRasterBandStats.statsGathered = QgsRasterBandStats::All;
244 : 0 : mStatistics.append( myRasterBandStats );
245 : :
246 : 0 : return myRasterBandStats;
247 : 0 : }
248 : :
249 : 0 : void QgsRasterInterface::initHistogram( QgsRasterHistogram &histogram,
250 : : int bandNo,
251 : : int binCount,
252 : : double minimum, double maximum,
253 : : const QgsRectangle &boundingBox,
254 : : int sampleSize,
255 : : bool includeOutOfRange )
256 : : {
257 : 0 : histogram.bandNumber = bandNo;
258 : 0 : histogram.minimum = minimum;
259 : 0 : histogram.maximum = maximum;
260 : 0 : histogram.includeOutOfRange = includeOutOfRange;
261 : :
262 : 0 : int mySrcDataType = sourceDataType( bandNo );
263 : :
264 : 0 : if ( std::isnan( histogram.minimum ) )
265 : : {
266 : : // TODO: this was OK when stats/histogram were calced in provider,
267 : : // but what TODO in other interfaces? Check for mInput for now.
268 : 0 : if ( !mInput && mySrcDataType == Qgis::Byte )
269 : : {
270 : 0 : histogram.minimum = 0; // see histogram() for shift for rounding
271 : 0 : }
272 : : else
273 : : {
274 : : // We need statistics -> avoid histogramDefaults in hasHistogram if possible
275 : : // TODO: use approximated statistics if approximated histogram is requested
276 : : // (theSampleSize > 0)
277 : 0 : QgsRasterBandStats stats = bandStatistics( bandNo, QgsRasterBandStats::Min, boundingBox, sampleSize );
278 : 0 : histogram.minimum = stats.minimumValue;
279 : : }
280 : 0 : }
281 : 0 : if ( std::isnan( histogram.maximum ) )
282 : : {
283 : 0 : if ( !mInput && mySrcDataType == Qgis::Byte )
284 : : {
285 : 0 : histogram.maximum = 255;
286 : 0 : }
287 : : else
288 : : {
289 : 0 : QgsRasterBandStats stats = bandStatistics( bandNo, QgsRasterBandStats::Max, boundingBox, sampleSize );
290 : 0 : histogram.maximum = stats.maximumValue;
291 : : }
292 : 0 : }
293 : :
294 : 0 : QgsRectangle finalExtent;
295 : 0 : if ( boundingBox.isEmpty() )
296 : : {
297 : 0 : finalExtent = extent();
298 : 0 : }
299 : : else
300 : : {
301 : 0 : finalExtent = extent().intersect( boundingBox );
302 : : }
303 : 0 : histogram.extent = finalExtent;
304 : :
305 : 0 : if ( sampleSize > 0 )
306 : : {
307 : : // Calc resolution from theSampleSize
308 : : double xRes, yRes;
309 : 0 : xRes = yRes = std::sqrt( ( static_cast<double>( finalExtent.width( ) ) * finalExtent.height() ) / sampleSize );
310 : :
311 : : // But limit by physical resolution
312 : 0 : if ( capabilities() & Size )
313 : : {
314 : 0 : double srcXRes = extent().width() / xSize();
315 : 0 : double srcYRes = extent().height() / ySize();
316 : 0 : if ( xRes < srcXRes ) xRes = srcXRes;
317 : 0 : if ( yRes < srcYRes ) yRes = srcYRes;
318 : 0 : }
319 : 0 : QgsDebugMsgLevel( QStringLiteral( "xRes = %1 yRes = %2" ).arg( xRes ).arg( yRes ), 4 );
320 : :
321 : 0 : histogram.width = static_cast <int>( finalExtent.width() / xRes );
322 : 0 : histogram.height = static_cast <int>( finalExtent.height() / yRes );
323 : 0 : }
324 : : else
325 : : {
326 : 0 : if ( capabilities() & Size )
327 : : {
328 : 0 : histogram.width = xSize();
329 : 0 : histogram.height = ySize();
330 : 0 : }
331 : : else
332 : : {
333 : 0 : histogram.width = 1000;
334 : 0 : histogram.height = 1000;
335 : : }
336 : : }
337 : 0 : QgsDebugMsgLevel( QStringLiteral( "theHistogram.width = %1 histogram.height = %2" ).arg( histogram.width ).arg( histogram.height ), 4 );
338 : :
339 : 0 : qint64 myBinCount = binCount;
340 : 0 : if ( myBinCount == 0 )
341 : : {
342 : : // TODO: this was OK when stats/histogram were calced in provider,
343 : : // but what TODO in other interfaces? Check for mInput for now.
344 : 0 : if ( !mInput && mySrcDataType == Qgis::Byte )
345 : : {
346 : 0 : myBinCount = 256; // Cannot store more values in byte
347 : 0 : }
348 : : else
349 : : {
350 : : // There is no best default value, to display something reasonable in histogram chart,
351 : : // binCount should be small, OTOH, to get precise data for cumulative cut, the number should be big.
352 : : // Because it is easier to define fixed lower value for the chart, we calc optimum binCount
353 : : // for higher resolution (to avoid calculating that where histogram() is used. In any case,
354 : : // it does not make sense to use more than width*height;
355 : :
356 : : // for Int16/Int32 make sure bin count <= actual range, because there is no sense in having
357 : : // bins at fractional values
358 : 0 : if ( !mInput && (
359 : 0 : mySrcDataType == Qgis::Int16 || mySrcDataType == Qgis::Int32 ||
360 : 0 : mySrcDataType == Qgis::UInt16 || mySrcDataType == Qgis::UInt32 ) )
361 : : {
362 : 0 : myBinCount = std::min( static_cast<qint64>( histogram.width ) * histogram.height, static_cast<qint64>( std::ceil( histogram.maximum - histogram.minimum + 1 ) ) );
363 : 0 : }
364 : : else
365 : : {
366 : : // This is for not integer types
367 : 0 : myBinCount = static_cast<qint64>( histogram.width ) * static_cast<qint64>( histogram.height );
368 : : }
369 : : }
370 : 0 : }
371 : : // Hard limit 10'000'000
372 : 0 : histogram.binCount = static_cast<int>( std::min( 10000000LL, myBinCount ) );
373 : 0 : QgsDebugMsgLevel( QStringLiteral( "theHistogram.binCount = %1" ).arg( histogram.binCount ), 4 );
374 : 0 : }
375 : :
376 : :
377 : 0 : bool QgsRasterInterface::hasHistogram( int bandNo,
378 : : int binCount,
379 : : double minimum, double maximum,
380 : : const QgsRectangle &extent,
381 : : int sampleSize,
382 : : bool includeOutOfRange )
383 : : {
384 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 binCount = %2 minimum = %3 maximum = %4 sampleSize = %5" ).arg( bandNo ).arg( binCount ).arg( minimum ).arg( maximum ).arg( sampleSize ), 4 );
385 : : // histogramDefaults() needs statistics if minimum or maximum is NaN ->
386 : : // do other checks which don't need statistics before histogramDefaults()
387 : 0 : if ( mHistograms.isEmpty() ) return false;
388 : :
389 : 0 : QgsRasterHistogram myHistogram;
390 : 0 : initHistogram( myHistogram, bandNo, binCount, minimum, maximum, extent, sampleSize, includeOutOfRange );
391 : :
392 : 0 : const auto constMHistograms = mHistograms;
393 : 0 : for ( const QgsRasterHistogram &histogram : constMHistograms )
394 : : {
395 : 0 : if ( histogram == myHistogram )
396 : : {
397 : 0 : QgsDebugMsgLevel( QStringLiteral( "Has cached histogram." ), 4 );
398 : 0 : return true;
399 : : }
400 : : }
401 : 0 : return false;
402 : 0 : }
403 : :
404 : 0 : QgsRasterHistogram QgsRasterInterface::histogram( int bandNo,
405 : : int binCount,
406 : : double minimum, double maximum,
407 : : const QgsRectangle &extent,
408 : : int sampleSize,
409 : : bool includeOutOfRange, QgsRasterBlockFeedback *feedback )
410 : : {
411 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 binCount = %2 minimum = %3 maximum = %4 sampleSize = %5" ).arg( bandNo ).arg( binCount ).arg( minimum ).arg( maximum ).arg( sampleSize ), 4 );
412 : :
413 : 0 : QgsRasterHistogram myHistogram;
414 : 0 : initHistogram( myHistogram, bandNo, binCount, minimum, maximum, extent, sampleSize, includeOutOfRange );
415 : :
416 : : // Find cached
417 : 0 : const auto constMHistograms = mHistograms;
418 : 0 : for ( const QgsRasterHistogram &histogram : constMHistograms )
419 : : {
420 : 0 : if ( histogram == myHistogram )
421 : : {
422 : 0 : QgsDebugMsgLevel( QStringLiteral( "Using cached histogram." ), 4 );
423 : 0 : return histogram;
424 : : }
425 : : }
426 : :
427 : 0 : int myBinCount = myHistogram.binCount;
428 : 0 : int myWidth = myHistogram.width;
429 : 0 : int myHeight = myHistogram.height;
430 : 0 : QgsRectangle myExtent = myHistogram.extent;
431 : 0 : myHistogram.histogramVector.resize( myBinCount );
432 : :
433 : 0 : int myXBlockSize = xBlockSize();
434 : 0 : int myYBlockSize = yBlockSize();
435 : 0 : if ( myXBlockSize == 0 ) // should not happen, but happens
436 : : {
437 : 0 : myXBlockSize = 500;
438 : 0 : }
439 : 0 : if ( myYBlockSize == 0 ) // should not happen, but happens
440 : : {
441 : 0 : myYBlockSize = 500;
442 : 0 : }
443 : :
444 : 0 : int myNXBlocks = ( myWidth + myXBlockSize - 1 ) / myXBlockSize;
445 : 0 : int myNYBlocks = ( myHeight + myYBlockSize - 1 ) / myYBlockSize;
446 : :
447 : 0 : double myXRes = myExtent.width() / myWidth;
448 : 0 : double myYRes = myExtent.height() / myHeight;
449 : :
450 : 0 : double myMinimum = myHistogram.minimum;
451 : 0 : double myMaximum = myHistogram.maximum;
452 : :
453 : : // To avoid rounding errors
454 : : // TODO: check this
455 : 0 : double myerval = ( myMaximum - myMinimum ) / myHistogram.binCount;
456 : 0 : myMinimum -= 0.1 * myerval;
457 : 0 : myMaximum += 0.1 * myerval;
458 : :
459 : 0 : QgsDebugMsgLevel( QStringLiteral( "binCount = %1 myMinimum = %2 myMaximum = %3" ).arg( myHistogram.binCount ).arg( myMinimum ).arg( myMaximum ), 4 );
460 : :
461 : 0 : double myBinSize = ( myMaximum - myMinimum ) / myBinCount;
462 : :
463 : : // TODO: progress signals
464 : 0 : bool isNoData = false;
465 : 0 : for ( int myYBlock = 0; myYBlock < myNYBlocks; myYBlock++ )
466 : : {
467 : 0 : for ( int myXBlock = 0; myXBlock < myNXBlocks; myXBlock++ )
468 : : {
469 : 0 : if ( feedback && feedback->isCanceled() )
470 : 0 : return myHistogram;
471 : :
472 : 0 : int myBlockWidth = std::min( myXBlockSize, myWidth - myXBlock * myXBlockSize );
473 : 0 : int myBlockHeight = std::min( myYBlockSize, myHeight - myYBlock * myYBlockSize );
474 : :
475 : 0 : double xmin = myExtent.xMinimum() + myXBlock * myXBlockSize * myXRes;
476 : 0 : double xmax = xmin + myBlockWidth * myXRes;
477 : 0 : double ymin = myExtent.yMaximum() - myYBlock * myYBlockSize * myYRes;
478 : 0 : double ymax = ymin - myBlockHeight * myYRes;
479 : :
480 : 0 : QgsRectangle myPartExtent( xmin, ymin, xmax, ymax );
481 : :
482 : 0 : std::unique_ptr< QgsRasterBlock > blk( block( bandNo, myPartExtent, myBlockWidth, myBlockHeight, feedback ) );
483 : :
484 : : // Collect the histogram counts.
485 : 0 : for ( qgssize i = 0; i < ( static_cast< qgssize >( myBlockHeight ) ) * myBlockWidth; i++ )
486 : : {
487 : 0 : double myValue = blk->valueAndNoData( i, isNoData );
488 : 0 : if ( isNoData )
489 : : {
490 : 0 : continue; // NULL
491 : : }
492 : :
493 : 0 : int myBinIndex = static_cast <int>( std::floor( ( myValue - myMinimum ) / myBinSize ) );
494 : :
495 : 0 : if ( ( myBinIndex < 0 || myBinIndex > ( myBinCount - 1 ) ) && !includeOutOfRange )
496 : : {
497 : 0 : continue;
498 : : }
499 : 0 : if ( myBinIndex < 0 ) myBinIndex = 0;
500 : 0 : if ( myBinIndex > ( myBinCount - 1 ) ) myBinIndex = myBinCount - 1;
501 : :
502 : 0 : myHistogram.histogramVector[myBinIndex] += 1;
503 : 0 : myHistogram.nonNullCount++;
504 : 0 : }
505 : 0 : }
506 : 0 : }
507 : :
508 : 0 : myHistogram.valid = true;
509 : 0 : mHistograms.append( myHistogram );
510 : :
511 : : #ifdef QGISDEBUG
512 : : QString hist;
513 : : for ( int i = 0; i < std::min( myHistogram.histogramVector.size(), 500 ); i++ )
514 : : {
515 : : hist += QString::number( myHistogram.histogramVector.value( i ) ) + ' ';
516 : : }
517 : : QgsDebugMsgLevel( QStringLiteral( "Histogram (max first 500 bins): " ) + hist, 4 );
518 : : #endif
519 : :
520 : 0 : return myHistogram;
521 : 0 : }
522 : :
523 : 0 : void QgsRasterInterface::cumulativeCut( int bandNo,
524 : : double lowerCount, double upperCount,
525 : : double &lowerValue, double &upperValue,
526 : : const QgsRectangle &extent,
527 : : int sampleSize )
528 : : {
529 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 lowerCount = %2 upperCount = %3 sampleSize = %4" ).arg( bandNo ).arg( lowerCount ).arg( upperCount ).arg( sampleSize ), 4 );
530 : :
531 : 0 : int mySrcDataType = sourceDataType( bandNo );
532 : :
533 : : // Init to NaN is better than histogram min/max to catch errors
534 : 0 : lowerValue = std::numeric_limits<double>::quiet_NaN();
535 : 0 : upperValue = std::numeric_limits<double>::quiet_NaN();
536 : :
537 : : //get band stats to specify real histogram min/max (fix #9793 Byte bands)
538 : 0 : QgsRasterBandStats stats = bandStatistics( bandNo, QgsRasterBandStats::Min, extent, sampleSize );
539 : 0 : if ( stats.maximumValue < stats.minimumValue )
540 : 0 : return;
541 : :
542 : : // for byte bands make sure bin count == actual range
543 : 0 : int myBinCount = ( mySrcDataType == Qgis::Byte ) ? int( std::ceil( stats.maximumValue - stats.minimumValue + 1 ) ) : 0;
544 : 0 : QgsRasterHistogram myHistogram = histogram( bandNo, myBinCount, stats.minimumValue, stats.maximumValue, extent, sampleSize );
545 : : //QgsRasterHistogram myHistogram = histogram( bandNo, 0, std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), extent, sampleSize );
546 : :
547 : 0 : double myBinXStep = ( myHistogram.maximum - myHistogram.minimum ) / myHistogram.binCount;
548 : 0 : int myCount = 0;
549 : 0 : int myMinCount = static_cast< int >( std::round( lowerCount * myHistogram.nonNullCount ) );
550 : 0 : int myMaxCount = static_cast< int >( std::round( upperCount * myHistogram.nonNullCount ) );
551 : 0 : bool myLowerFound = false;
552 : 0 : QgsDebugMsgLevel( QStringLiteral( "binCount = %1 minimum = %2 maximum = %3 myBinXStep = %4" ).arg( myHistogram.binCount ).arg( myHistogram.minimum ).arg( myHistogram.maximum ).arg( myBinXStep ), 4 );
553 : 0 : QgsDebugMsgLevel( QStringLiteral( "myMinCount = %1 myMaxCount = %2" ).arg( myMinCount ).arg( myMaxCount ), 4 );
554 : :
555 : 0 : for ( int myBin = 0; myBin < myHistogram.histogramVector.size(); myBin++ )
556 : : {
557 : 0 : int myBinValue = myHistogram.histogramVector.value( myBin );
558 : 0 : myCount += myBinValue;
559 : 0 : if ( !myLowerFound && myCount > myMinCount )
560 : : {
561 : 0 : lowerValue = myHistogram.minimum + myBin * myBinXStep;
562 : 0 : myLowerFound = true;
563 : 0 : QgsDebugMsgLevel( QStringLiteral( "found lowerValue %1 at bin %2" ).arg( lowerValue ).arg( myBin ), 4 );
564 : 0 : }
565 : 0 : if ( myCount >= myMaxCount )
566 : : {
567 : 0 : upperValue = myHistogram.minimum + myBin * myBinXStep;
568 : 0 : QgsDebugMsgLevel( QStringLiteral( "found upperValue %1 at bin %2" ).arg( upperValue ).arg( myBin ), 4 );
569 : 0 : break;
570 : : }
571 : 0 : }
572 : :
573 : : // fix integer data - round down/up
574 : 0 : if ( mySrcDataType == Qgis::Byte ||
575 : 0 : mySrcDataType == Qgis::Int16 || mySrcDataType == Qgis::Int32 ||
576 : 0 : mySrcDataType == Qgis::UInt16 || mySrcDataType == Qgis::UInt32 )
577 : : {
578 : 0 : if ( !std::isnan( lowerValue ) )
579 : 0 : lowerValue = std::floor( lowerValue );
580 : 0 : if ( !std::isnan( upperValue ) )
581 : 0 : upperValue = std::ceil( upperValue );
582 : 0 : }
583 : 0 : }
584 : :
585 : 0 : QString QgsRasterInterface::capabilitiesString() const
586 : : {
587 : 0 : QStringList abilitiesList;
588 : :
589 : 0 : int abilities = capabilities();
590 : :
591 : : // Not all all capabilities are here (Size, IdentifyValue, IdentifyText,
592 : : // IdentifyHtml, IdentifyFeature) because those are quite technical and probably
593 : : // would be confusing for users
594 : :
595 : 0 : if ( abilities & QgsRasterInterface::Identify )
596 : : {
597 : 0 : abilitiesList += tr( "Identify" );
598 : 0 : }
599 : :
600 : 0 : if ( abilities & QgsRasterInterface::Create )
601 : : {
602 : 0 : abilitiesList += tr( "Create Datasources" );
603 : 0 : }
604 : :
605 : 0 : if ( abilities & QgsRasterInterface::Remove )
606 : : {
607 : 0 : abilitiesList += tr( "Remove Datasources" );
608 : 0 : }
609 : :
610 : 0 : if ( abilities & QgsRasterInterface::BuildPyramids )
611 : : {
612 : 0 : abilitiesList += tr( "Build Pyramids" );
613 : 0 : }
614 : :
615 : 0 : QgsDebugMsgLevel( "Capability: " + abilitiesList.join( QLatin1String( ", " ) ), 4 );
616 : :
617 : 0 : return abilitiesList.join( QLatin1String( ", " ) );
618 : 0 : }
619 : :
620 : 0 : QString QgsRasterInterface::generateBandName( int bandNumber ) const
621 : : {
622 : 0 : if ( mInput )
623 : 0 : return mInput->generateBandName( bandNumber );
624 : :
625 : 0 : return tr( "Band" ) + QStringLiteral( " %1" ) .arg( bandNumber, 1 + static_cast< int >( std::log10( static_cast< double >( bandCount() ) ) ), 10, QChar( '0' ) );
626 : 0 : }
627 : :
628 : 0 : QString QgsRasterInterface::colorInterpretationName( int bandNo ) const
629 : : {
630 : 0 : if ( mInput )
631 : 0 : return mInput->colorInterpretationName( bandNo );
632 : :
633 : 0 : return QString();
634 : 0 : }
635 : :
636 : 0 : QString QgsRasterInterface::displayBandName( int bandNumber ) const
637 : : {
638 : 0 : QString name = generateBandName( bandNumber );
639 : 0 : QString colorInterp = colorInterpretationName( bandNumber );
640 : 0 : if ( colorInterp != QLatin1String( "Undefined" ) )
641 : : {
642 : 0 : name.append( QStringLiteral( " (%1)" ).arg( colorInterp ) );
643 : 0 : }
644 : 0 : return name;
645 : 0 : }
|
