From 4db175e29abd61beb3db7b3504569381f8417ff3 Mon Sep 17 00:00:00 2001 From: Tim Sutton Date: Thu, 8 Sep 2011 23:04:52 +0200 Subject: [PATCH] Backported raster speed up (using providers to compute stats natively) to release branch. --- src/core/qgsrasterdataprovider.cpp | 203 +++++++++++++++++++++-- src/core/qgsrasterdataprovider.h | 41 +++-- src/providers/gdal/qgsgdalprovider.cpp | 213 +++++++++++++++++-------- src/providers/gdal/qgsgdalprovider.h | 23 ++- src/providers/wms/CMakeLists.txt | 6 +- 5 files changed, 385 insertions(+), 101 deletions(-) diff --git a/src/core/qgsrasterdataprovider.cpp b/src/core/qgsrasterdataprovider.cpp index b707fe6df307..b7206ce4474c 100644 --- a/src/core/qgsrasterdataprovider.cpp +++ b/src/core/qgsrasterdataprovider.cpp @@ -14,7 +14,6 @@ * (at your option) any later version. * * * ***************************************************************************/ -/* $Id$ */ #include "qgsrasterdataprovider.h" #include "qgsrasterprojector.h" @@ -54,7 +53,8 @@ void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle const & viewExt // TODO: init data by nulls // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers - if ( myProjector.srcRows() <= 0 || myProjector.srcCols() <= 0 ) return; + if ( myProjector.srcRows() <= 0 || myProjector.srcCols() <= 0 ) + return; // Allocate memory for not projected source data int mySize = dataTypeSize( bandNo ) / 8; @@ -174,6 +174,7 @@ QString QgsRasterDataProvider::capabilitiesString() const bool QgsRasterDataProvider::identify( const QgsPoint& thePoint, QMap& theResults ) { + Q_UNUSED( thePoint ); theResults.clear(); return false; } @@ -200,31 +201,31 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo ) double d; switch ( type ) { - case QgsRasterDataProvider::Byte: + case Byte: uc = ( unsigned char )noval; memcpy( data, &uc, size ); break; - case QgsRasterDataProvider::UInt16: + case UInt16: us = ( unsigned short )noval; memcpy( data, &us, size ); break; - case QgsRasterDataProvider::Int16: + case Int16: s = ( short )noval; memcpy( data, &s, size ); break; - case QgsRasterDataProvider::UInt32: + case UInt32: ui = ( unsigned int )noval; memcpy( data, &ui, size ); break; - case QgsRasterDataProvider::Int32: + case Int32: i = ( int )noval; memcpy( data, &i, size ); break; - case QgsRasterDataProvider::Float32: + case Float32: f = ( float )noval; memcpy( data, &f, size ); break; - case QgsRasterDataProvider::Float64: + case Float64: d = ( double )noval; memcpy( data, &d, size ); break; @@ -234,4 +235,188 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo ) return ba; } + +QgsRasterBandStats QgsRasterDataProvider::bandStatistics( int theBandNo ) +{ + double myNoDataValue = noDataValue(); + QgsRasterBandStats myRasterBandStats; + myRasterBandStats.elementCount = 0; // because we'll be counting only VALID pixels later + + int myDataType = dataType( theBandNo ); + + int myNXBlocks, myNYBlocks, myXBlockSize, myYBlockSize; + myXBlockSize = xBlockSize(); + myYBlockSize = yBlockSize(); + + myNXBlocks = ( xSize() + myXBlockSize - 1 ) / myXBlockSize; + myNYBlocks = ( ySize() + myYBlockSize - 1 ) / myYBlockSize; + + void *myData = CPLMalloc( myXBlockSize * myYBlockSize * ( dataTypeSize( theBandNo ) / 8 ) ); + + // unfortunately we need to make two passes through the data to calculate stddev + bool myFirstIterationFlag = true; + + int myBandXSize = xSize(); + int myBandYSize = ySize(); + for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ ) + { + for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ ) + { + int nXValid, nYValid; + readBlock( theBandNo, iXBlock, iYBlock, myData ); + + // Compute the portion of the block that is valid + // for partial edge blocks. + if (( iXBlock + 1 ) * myXBlockSize > myBandXSize ) + nXValid = myBandXSize - iXBlock * myXBlockSize; + else + nXValid = myXBlockSize; + + if (( iYBlock + 1 ) * myYBlockSize > myBandYSize ) + nYValid = myBandYSize - iYBlock * myYBlockSize; + else + nYValid = myYBlockSize; + + // Collect the histogram counts. + for ( int iY = 0; iY < nYValid; iY++ ) + { + for ( int iX = 0; iX < nXValid; iX++ ) + { + double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) ); + //QgsDebugMsg ( QString ( "%1 %2 value %3" ).arg (iX).arg(iY).arg( myValue ) ); + + if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) ) + { + continue; // NULL + } + + myRasterBandStats.sum += myValue; + ++myRasterBandStats.elementCount; + //only use this element if we have a non null element + if ( myFirstIterationFlag ) + { + //this is the first iteration so initialise vars + myFirstIterationFlag = false; + myRasterBandStats.minimumValue = myValue; + myRasterBandStats.maximumValue = myValue; + } //end of true part for first iteration check + else + { + //this is done for all subsequent iterations + if ( myValue < myRasterBandStats.minimumValue ) + { + myRasterBandStats.minimumValue = myValue; + } + if ( myValue > myRasterBandStats.maximumValue ) + { + myRasterBandStats.maximumValue = myValue; + } + } //end of false part for first iteration check + } + } + } //end of column wise loop + } //end of row wise loop + + + //end of first pass through data now calculate the range + myRasterBandStats.range = myRasterBandStats.maximumValue - myRasterBandStats.minimumValue; + //calculate the mean + myRasterBandStats.mean = myRasterBandStats.sum / myRasterBandStats.elementCount; + + //for the second pass we will get the sum of the squares / mean + for ( int iYBlock = 0; iYBlock < myNYBlocks; iYBlock++ ) + { + for ( int iXBlock = 0; iXBlock < myNXBlocks; iXBlock++ ) + { + int nXValid, nYValid; + + readBlock( theBandNo, iXBlock, iYBlock, myData ); + + // Compute the portion of the block that is valid + // for partial edge blocks. + if (( iXBlock + 1 ) * myXBlockSize > myBandXSize ) + nXValid = myBandXSize - iXBlock * myXBlockSize; + else + nXValid = myXBlockSize; + + if (( iYBlock + 1 ) * myYBlockSize > myBandYSize ) + nYValid = myBandYSize - iYBlock * myYBlockSize; + else + nYValid = myYBlockSize; + + // Collect the histogram counts. + for ( int iY = 0; iY < nYValid; iY++ ) + { + for ( int iX = 0; iX < nXValid; iX++ ) + { + double myValue = readValue( myData, myDataType, iX + ( iY * myXBlockSize ) ); + //QgsDebugMsg ( "myValue = " + QString::number(myValue) ); + + if ( mValidNoDataValue && ( qAbs( myValue - myNoDataValue ) <= TINY_VALUE ) ) + { + continue; // NULL + } + + myRasterBandStats.sumOfSquares += static_cast < double > + ( pow( myValue - myRasterBandStats.mean, 2 ) ); + } + } + } //end of column wise loop + } //end of row wise loop + + //divide result by sample size - 1 and get square root to get stdev + myRasterBandStats.stdDev = static_cast < double >( sqrt( myRasterBandStats.sumOfSquares / + ( myRasterBandStats.elementCount - 1 ) ) ); + +#ifdef QGISDEBUG + QgsLogger::debug( "************ STATS **************", 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "VALID NODATA", mValidNoDataValue, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "NULL", noDataValue() , 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "MIN", myRasterBandStats.minimumValue, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "MAX", myRasterBandStats.maximumValue, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "RANGE", myRasterBandStats.range, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "MEAN", myRasterBandStats.mean, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "STDDEV", myRasterBandStats.stdDev, 1, __FILE__, __FUNCTION__, __LINE__ ); +#endif + + CPLFree( myData ); + myRasterBandStats.statsGathered = true; + return myRasterBandStats; +} + +double QgsRasterDataProvider::readValue( void *data, int type, int index ) +{ + if ( !data ) + return mValidNoDataValue ? noDataValue() : 0.0; + + switch ( type ) + { + case Byte: + return ( double )(( GByte * )data )[index]; + break; + case UInt16: + return ( double )(( GUInt16 * )data )[index]; + break; + case Int16: + return ( double )(( GInt16 * )data )[index]; + break; + case UInt32: + return ( double )(( GUInt32 * )data )[index]; + break; + case Int32: + return ( double )(( GInt32 * )data )[index]; + break; + case Float32: + return ( double )(( float * )data )[index]; + break; + case Float64: + return ( double )(( double * )data )[index]; + break; + default: + QgsLogger::warning( "GDAL data type is not supported" ); + } + + return mValidNoDataValue ? noDataValue() : 0.0; +} + // ENDS diff --git a/src/core/qgsrasterdataprovider.h b/src/core/qgsrasterdataprovider.h index a72701ad97dd..5d54b9869486 100644 --- a/src/core/qgsrasterdataprovider.h +++ b/src/core/qgsrasterdataprovider.h @@ -14,7 +14,6 @@ * (at your option) any later version. * * * ***************************************************************************/ -/* $Id$ */ /* Thank you to Marco Hugentobler for the original vector DataProvider */ @@ -29,12 +28,12 @@ #include "qgscolorrampshader.h" #include "qgsrasterpyramid.h" #include "qgscoordinatereferencesystem.h" - +#include "qgsrasterbandstats.h" +#include "cpl_conv.h" #include class QImage; class QgsPoint; -class QgsRasterBandStats; class QByteArray; #define TINY_VALUE std::numeric_limits::epsilon() * 20 @@ -110,10 +109,11 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider }; // Progress types - enum Progress + enum RasterProgressType { ProgressHistogram = 0, - ProgressPyramids = 1 + ProgressPyramids = 1, + ProgressStatistics = 2 }; QgsRasterDataProvider(); @@ -184,6 +184,7 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider */ virtual int srcDataType( int bandNo ) const { + Q_UNUSED( bandNo ); return QgsRasterDataProvider::UnknownDataType; } @@ -234,6 +235,7 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider /** Returns data type for the band specified by number */ virtual int colorInterpretation( int theBandNo ) const { + Q_UNUSED( theBandNo ); return QgsRasterDataProvider::UndefinedColorInterpretation; } @@ -315,21 +317,27 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider /** read block of data */ // TODO clarify what happens on the last block (the part outside raster) - virtual void readBlock( int bandNo, int xBlock, int yBlock, void *data ) {} + virtual void readBlock( int bandNo, int xBlock, int yBlock, void *data ) + { Q_UNUSED( bandNo ); Q_UNUSED( xBlock ); Q_UNUSED( yBlock ); Q_UNUSED( data ); } /** read block of data using give extent and size */ - virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, void *data ) {}; + virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, void *data ) + { Q_UNUSED( bandNo ); Q_UNUSED( viewExtent ); Q_UNUSED( width ); Q_UNUSED( height ); Q_UNUSED( data ); } /** read block of data using give extent and size */ virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, QgsCoordinateReferenceSystem theSrcCRS, QgsCoordinateReferenceSystem theDestCRS, void *data ); + + /* Read a value from a data block at a given index. */ + virtual double readValue( void *data, int type, int index ); /** value representing null data */ virtual double noDataValue() const { return 0; } - virtual double minimumValue( int bandNo )const { return 0; } - virtual double maximumValue( int bandNo )const { return 0; } + virtual double minimumValue( int bandNo ) const { Q_UNUSED( bandNo ); return 0; } + virtual double maximumValue( int bandNo ) const { Q_UNUSED( bandNo ); return 0; } - virtual QList colorTable( int bandNo )const { return QList(); } + virtual QList colorTable( int bandNo ) const + { Q_UNUSED( bandNo ); return QList(); } // Defined in parent /** \brief Returns the sublayers of this layer - Useful for providers that manage their own layers, such as WMS */ @@ -345,12 +353,14 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider int theBinCountInt = 256, bool theIgnoreOutOfRangeFlag = true, bool theThoroughBandScanFlag = false - ) {}; + ) + { Q_UNUSED( theBandNoInt ); Q_UNUSED( theBandStats ); Q_UNUSED( theBinCountInt ); Q_UNUSED( theIgnoreOutOfRangeFlag ); Q_UNUSED( theThoroughBandScanFlag ); } /** \brief Create pyramid overviews */ virtual QString buildPyramids( const QList & thePyramidList, const QString & theResamplingMethod = "NEAREST", - bool theTryInternalFlag = false ) { return "FAILED_NOT_SUPPORTED"; }; + bool theTryInternalFlag = false ) + { Q_UNUSED( thePyramidList ); Q_UNUSED( theResamplingMethod ); Q_UNUSED( theTryInternalFlag ); return "FAILED_NOT_SUPPORTED"; }; /** \brief Accessor for ths raster layers pyramid list. A pyramid list defines the * POTENTIAL pyramids that can be in a raster. To know which of the pyramid layers @@ -359,13 +369,16 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider */ virtual QList buildPyramidList() { return QList(); }; + /** If the provider supports it, return band stats for the + given band. Default behaviour is to blockwise read the data + and generate the stats unless the provider overloads this function. */ + virtual QgsRasterBandStats bandStatistics( int theBandNo ); /** \brief helper function to create zero padded band names */ QString generateBandName( int theBandNumber ) { return tr( "Band" ) + QString( " %1" ) .arg( theBandNumber, 1 + ( int ) log10(( float ) bandCount() ), 10, QChar( '0' ) ); - } - + }; /** * Get metadata in a format suitable for feeding directly diff --git a/src/providers/gdal/qgsgdalprovider.cpp b/src/providers/gdal/qgsgdalprovider.cpp index 57d1851c42db..3cba3ba2f788 100644 --- a/src/providers/gdal/qgsgdalprovider.cpp +++ b/src/providers/gdal/qgsgdalprovider.cpp @@ -16,7 +16,6 @@ * * ***************************************************************************/ -/* $Id: qgsgdalprovider.cpp 11522 2009-08-28 14:49:22Z jef $ */ #include "qgslogger.h" #include "qgsgdalprovider.h" @@ -27,6 +26,7 @@ #include "qgsrectangle.h" #include "qgscoordinatereferencesystem.h" #include "qgsrasterbandstats.h" +#include "qgsrasterlayer.h" #include "qgsrasterpyramid.h" #include @@ -103,7 +103,10 @@ QgsGdalProvider::QgsGdalProvider( QString const & uri ) registerGdalDrivers(); // To get buildSupportedRasterFileFilter the provider is called with empty uri - if ( uri.isEmpty() ) return; + if ( uri.isEmpty() ) + { + return; + } mGdalDataset = NULL; @@ -285,51 +288,6 @@ QgsGdalProvider::QgsGdalProvider( QString const & uri ) QgsDebugMsg( QString( "mNoDataValue[%1] = %1" ).arg( i - 1 ).arg( mNoDataValue[i-1] ) ); } - // This block of code was in old version in QgsRasterLayer::bandStatistics - //ifdefs below to remove compiler warning about unused vars -#ifdef QGISDEBUG -#if 0 - int success; - double GDALminimum = GDALGetRasterMinimum( myGdalBand, &success ); - - if ( ! success ) - { - QgsDebugMsg( "myGdalBand->GetMinimum() failed" ); - } - - double GDALmaximum = GDALGetRasterMaximum( myGdalBand, &success ); - - if ( ! success ) - { - QgsDebugMsg( "myGdalBand->GetMaximum() failed" ); - } - - double GDALnodata = GDALGetRasterNoDataValue( myGdalBand, &success ); - - if ( ! success ) - { - QgsDebugMsg( "myGdalBand->GetNoDataValue() failed" ); - } - - QgsLogger::debug( "GDALminium: ", GDALminimum, __FILE__, __FUNCTION__, __LINE__ ); - QgsLogger::debug( "GDALmaximum: ", GDALmaximum, __FILE__, __FUNCTION__, __LINE__ ); - QgsLogger::debug( "GDALnodata: ", GDALnodata, __FILE__, __FUNCTION__, __LINE__ ); - - double GDALrange[2]; // calculated min/max, as opposed to the - // dataset provided - - GDALComputeRasterMinMax( myGdalBand, 1, GDALrange ); - QgsLogger::debug( "approximate computed GDALminium:", GDALrange[0], __FILE__, __FUNCTION__, __LINE__, 1 ); - QgsLogger::debug( "approximate computed GDALmaximum:", GDALrange[1], __FILE__, __FUNCTION__, __LINE__, 1 ); - - GDALComputeRasterMinMax( myGdalBand, 0, GDALrange ); - QgsLogger::debug( "exactly computed GDALminium:", GDALrange[0] ); - QgsLogger::debug( "exactly computed GDALmaximum:", GDALrange[1] ); - - QgsDebugMsg( "starting manual stat computation" ); -#endif -#endif - mValid = true; QgsDebugMsg( "end" ); } @@ -388,7 +346,10 @@ QgsGdalProvider::~QgsGdalProvider() // This was used by raster layer to reload data void QgsGdalProvider::closeDataset() { - if ( !mValid ) return; + if ( !mValid ) + { + return; + } mValid = false; GDALDereferenceDataset( mGdalBaseDataset ); @@ -520,6 +481,7 @@ QString QgsGdalProvider::metadata() // Not supported by GDAL QImage* QgsGdalProvider::draw( QgsRectangle const & viewExtent, int pixelWidth, int pixelHeight ) { + Q_UNUSED( viewExtent ); QgsDebugMsg( "pixelWidth = " + QString::number( pixelWidth ) ); QgsDebugMsg( "pixelHeight = " + QString::number( pixelHeight ) ); QgsDebugMsg( "viewExtent: " + viewExtent.toString() ); @@ -916,25 +878,31 @@ double QgsGdalProvider::noDataValue() const void QgsGdalProvider::computeMinMax( int theBandNo ) { QgsDebugMsg( QString( "theBandNo = %1 mMinMaxComputed = %2" ).arg( theBandNo ).arg( mMinMaxComputed[theBandNo-1] ) ); - if ( mMinMaxComputed[theBandNo-1] ) return; - double GDALrange[2]; + if ( mMinMaxComputed[theBandNo-1] ) + { + return; + } GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, theBandNo ); - GDALComputeRasterMinMax( myGdalBand, 1, GDALrange ); //Approximate - QgsDebugMsg( QString( "GDALrange[0] = %1 GDALrange[1] = %2" ).arg( GDALrange[0] ).arg( GDALrange[1] ) ); - mMinimum[theBandNo-1] = GDALrange[0]; - mMaximum[theBandNo-1] = GDALrange[1]; + int bGotMin, bGotMax; + double adfMinMax[2]; + adfMinMax[0] = GDALGetRasterMinimum( myGdalBand, &bGotMin ); + adfMinMax[1] = GDALGetRasterMaximum( myGdalBand, &bGotMax ); + if ( !( bGotMin && bGotMax ) ) + { + GDALComputeRasterMinMax( myGdalBand, TRUE, adfMinMax ); + } + mMinimum[theBandNo-1] = adfMinMax[0]; + mMaximum[theBandNo-1] = adfMinMax[1]; } double QgsGdalProvider::minimumValue( int theBandNo ) const { QgsDebugMsg( QString( "theBandNo = %1" ).arg( theBandNo ) ); - //computeMinMax ( theBandNo ); return mMinimum[theBandNo-1]; } double QgsGdalProvider::maximumValue( int theBandNo ) const { QgsDebugMsg( QString( "theBandNo = %1" ).arg( theBandNo ) ); - //computeMinMax ( theBandNo ); return mMaximum[theBandNo-1]; } @@ -991,8 +959,8 @@ QList QgsGdalProvider::colorTable( int theBan else if ( myColorInterpretation == GCI_PaletteIndex ) { QgsColorRampShader::ColorRampItem myColorRampItem; - myColorRampItem.label = ""; myColorRampItem.value = ( double )myIterator; + myColorRampItem.label = QString::number( myColorRampItem.value ); //Branch on palette interpretation if ( myPaletteInterpretation == GPI_RGB ) { @@ -1225,11 +1193,13 @@ bool QgsGdalProvider::isValid() QString QgsGdalProvider::identifyAsText( const QgsPoint& point ) { + Q_UNUSED( point ); return QString( "Not implemented" ); } QString QgsGdalProvider::identifyAsHtml( const QgsPoint& point ) { + Q_UNUSED( point ); return QString( "Not implemented" ); } @@ -1287,7 +1257,8 @@ void QgsGdalProvider::populateHistogram( int theBandNo, QgsRasterBandStats & t //vector is not equal to the number of bins //i.e if the histogram has never previously been generated or the user has //selected a new number of bins. - if ( theBandStats.histogramVector->size() != theBinCount || + if ( theBandStats.histogramVector == 0 || + theBandStats.histogramVector->size() != theBinCount || theIgnoreOutOfRangeFlag != theBandStats.isHistogramOutOfRange || theHistogramEstimatedFlag != theBandStats.isHistogramEstimated ) { @@ -1363,7 +1334,7 @@ QString QgsGdalProvider::buildPyramids( QList const & theRaste // TODO add signal and connect from rasterlayer //emit drawingProgress( 0, 0 ); - //first test if the file is writeable + //first test if the file is writable //QFileInfo myQFile( mDataSource ); QFileInfo myQFile( dataSourceUri() ); @@ -1633,17 +1604,7 @@ QGISEXTERN bool isProvider() return true; } -/** - Builds the list of file filter strings to later be used by - QgisApp::addRasterLayer() - - We query GDAL for a list of supported raster formats; we then build - a list of file filter strings from that list. We return a string - that contains this list that is suitable for use in a - QFileDialog::getOpenFileNames() call. - -*/ -QGISEXTERN void buildSupportedRasterFileFilter( QString & theFileFiltersString ) +void buildSupportedRasterFileFilterAndExtensions( QString & theFileFiltersString, QStringList & theExtensions, QStringList & theWildcards ) { QgsDebugMsg( "Entered" ); // first get the GDAL driver manager @@ -1739,6 +1700,7 @@ QGISEXTERN void buildSupportedRasterFileFilter( QString & theFileFiltersString ) { // XXX add check for SDTS; in that case we want (*CATD.DDF) QString glob = "*." + myGdalDriverExtension.replace( "/", " *." ); + theExtensions << myGdalDriverExtension.replace( "/", "" ).replace( "*", "" ).replace( ".", "" ); // Add only the first JP2 driver found to the filter list (it's the one GDAL uses) if ( myGdalDriverDescription == "JPEG2000" || myGdalDriverDescription.startsWith( "JP2" ) ) // JP2ECW, JP2KAK, JP2MrSID @@ -1748,14 +1710,17 @@ QGISEXTERN void buildSupportedRasterFileFilter( QString & theFileFiltersString ) jp2Driver = myGdalDriver; // first JP2 driver found glob += " *.j2k"; // add alternate extension + theExtensions << "j2k"; } else if ( myGdalDriverDescription == "GTiff" ) { glob += " *.tiff"; + theExtensions << "tiff"; } else if ( myGdalDriverDescription == "JPEG" ) { glob += " *.jpeg"; + theExtensions << "jpeg"; } theFileFiltersString += ";;[GDAL] " + myGdalDriverLongName + " (" + glob.toLower() + " " + glob.toUpper() + ")"; @@ -1785,6 +1750,7 @@ QGISEXTERN void buildSupportedRasterFileFilter( QString & theFileFiltersString ) { QString glob = "*.dem"; theFileFiltersString += ";;[GDAL] " + myGdalDriverLongName + " (" + glob.toLower() + " " + glob.toUpper() + ")"; + theExtensions << "dem"; } else if ( myGdalDriverDescription.startsWith( "DTED" ) ) { @@ -1793,12 +1759,26 @@ QGISEXTERN void buildSupportedRasterFileFilter( QString & theFileFiltersString ) glob += " *.dt1"; glob += " *.dt2"; theFileFiltersString += ";;[GDAL] " + myGdalDriverLongName + " (" + glob.toLower() + " " + glob.toUpper() + ")"; + theExtensions << "dt0" << "dt1" << "dt2"; } else if ( myGdalDriverDescription.startsWith( "MrSID" ) ) { // MrSID use "*.sid" QString glob = "*.sid"; theFileFiltersString += ";;[GDAL] " + myGdalDriverLongName + " (" + glob.toLower() + " " + glob.toUpper() + ")"; + theExtensions << "sid"; + } + else if ( myGdalDriverDescription.startsWith( "EHdr" ) ) + { + QString glob = "*.bil"; + theFileFiltersString += ";;[GDAL] " + myGdalDriverLongName + " (" + glob.toLower() + " " + glob.toUpper() + ")"; + theExtensions << "bil"; + } + else if ( myGdalDriverDescription.startsWith( "AIG" ) ) + { + QString glob = "hdr.adf"; + theFileFiltersString += ";;[GDAL] " + myGdalDriverLongName + " (" + glob.toLower() + " " + glob.toUpper() + ")"; + theWildcards << "hdr.adf"; } else { @@ -1848,3 +1828,96 @@ QGISEXTERN bool isValidRasterFileName( QString const & theFileNameQString, QStri return true; } } + + + +QgsRasterBandStats QgsGdalProvider::bandStatistics( int theBandNo ) +{ + GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, theBandNo ); + QgsRasterBandStats myRasterBandStats; + int bApproxOK = false; + double pdfMin; + double pdfMax; + double pdfMean; + double pdfStdDev; + QgsGdalProgress myProg; + myProg.type = ProgressHistogram; + myProg.provider = this; + + // double myerval = + // GDALComputeRasterStatistics ( + // myGdalBand, bApproxOK, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev, + // progressCallback, &myProg ) ; + // double myerval = + // GDALGetRasterStatistics ( myGdalBand, bApproxOK, TRUE, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev); + // double myerval = + // GDALGetRasterStatisticsProgress ( myGdalBand, bApproxOK, TRUE, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev, + // progressCallback, &myProg ); + + // try to fetch the cached stats (bForce=FALSE) + CPLErr myerval = + GDALGetRasterStatistics( myGdalBand, bApproxOK, FALSE, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev ); + + // if cached stats are not found, compute them + if ( CE_Warning == myerval ) + { + myerval = GDALComputeRasterStatistics( myGdalBand, bApproxOK, + &pdfMin, &pdfMax, &pdfMean, &pdfStdDev, + progressCallback, &myProg ) ; + } + + // if stats are found populate the QgsRasterBandStats object + if ( CE_None == myerval ) + { + + myRasterBandStats.bandName = generateBandName( theBandNo ); + myRasterBandStats.bandNumber = theBandNo; + myRasterBandStats.range = pdfMax - pdfMin; + myRasterBandStats.minimumValue = pdfMin; + myRasterBandStats.maximumValue = pdfMax; + //calculate the mean + myRasterBandStats.mean = pdfMean; + myRasterBandStats.sum = 0; //not available via gdal + myRasterBandStats.elementCount = mWidth * mHeight; + myRasterBandStats.sumOfSquares = 0; //not available via gdal + myRasterBandStats.stdDev = pdfStdDev; + myRasterBandStats.statsGathered = true; + +#ifdef QGISDEBUG + QgsLogger::debug( "************ STATS **************", 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "VALID NODATA", mValidNoDataValue, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "MIN", myRasterBandStats.minimumValue, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "MAX", myRasterBandStats.maximumValue, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "RANGE", myRasterBandStats.range, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "MEAN", myRasterBandStats.mean, 1, __FILE__, __FUNCTION__, __LINE__ ); + QgsLogger::debug( "STDDEV", myRasterBandStats.stdDev, 1, __FILE__, __FUNCTION__, __LINE__ ); +#endif + + myRasterBandStats.statsGathered = true; + + } + + return myRasterBandStats; + +} // QgsGdalProvider::bandStatistics + +/** + Builds the list of file filter strings to later be used by + QgisApp::addRasterLayer() + + We query GDAL for a list of supported raster formats; we then build + a list of file filter strings from that list. We return a string + that contains this list that is suitable for use in a + QFileDialog::getOpenFileNames() call. + +*/ +QGISEXTERN void buildSupportedRasterFileFilter( QString & theFileFiltersString ) +{ + QStringList exts; + QStringList wildcards; + buildSupportedRasterFileFilterAndExtensions( theFileFiltersString, exts, wildcards ); +} + + + + diff --git a/src/providers/gdal/qgsgdalprovider.h b/src/providers/gdal/qgsgdalprovider.h index 6992b220e31b..84fc9e038b3d 100644 --- a/src/providers/gdal/qgsgdalprovider.h +++ b/src/providers/gdal/qgsgdalprovider.h @@ -16,7 +16,6 @@ * * ***************************************************************************/ -/* $Id: qgsgdalprovider.h 12528 2009-12-20 12:29:07Z jef $ */ #ifndef QGSGDALPROVIDER_H #define QGSGDALPROVIDER_H @@ -26,6 +25,7 @@ #include "qgsrasterdataprovider.h" #include "qgsrectangle.h" #include "qgscolorrampshader.h" +#include "qgsrasterbandstats.h" #include #include @@ -208,17 +208,27 @@ class QgsGdalProvider : public QgsRasterDataProvider QString metadata(); // Following methods specific for WMS are not used at all in this provider and should be removed IMO from qgsdataprovider.h - void addLayers( QStringList const & layers, QStringList const & styles = QStringList() ) {} - QStringList supportedImageEncodings() { return QStringList();} + void addLayers( QStringList const &layers, QStringList const &styles = QStringList() ) + { Q_UNUSED( layers ); Q_UNUSED( styles ); } + QStringList supportedImageEncodings() { return QStringList(); } QString imageEncoding() const { return QString(); } - void setImageEncoding( QString const & mimeType ) {} - void setImageCrs( QString const & crs ) {} + void setImageEncoding( QString const &mimeType ) + { Q_UNUSED( mimeType ); } + void setImageCrs( QString const &crs ) + { Q_UNUSED( crs ); } /** \brief ensures that GDAL drivers are registered, but only once */ static void registerGdalDrivers(); /** \brief Returns the sublayers of this layer - Useful for providers that manage their own layers, such as WMS */ QStringList subLayers() const; + /** \brief If the provider supports it, return band stats for the + given band. + @note added in QGIS 1.7 + @note overloads virtual method from QgsRasterProvider::bandStatistics + + */ + QgsRasterBandStats bandStatistics( int theBandNo ); void populateHistogram( int theBandNoInt, QgsRasterBandStats & theBandStats, @@ -270,8 +280,6 @@ class QgsGdalProvider : public QgsRasterDataProvider // List of estimated max values, index 0 for band 1 QList mMaximum; - //GDALDataType mGdalDataType; - /** \brief Pointer to the gdaldataset */ GDALDatasetH mGdalBaseDataset; @@ -287,5 +295,6 @@ class QgsGdalProvider : public QgsRasterDataProvider }; + #endif diff --git a/src/providers/wms/CMakeLists.txt b/src/providers/wms/CMakeLists.txt index 1adc5d99f9d1..995aaefacf6b 100644 --- a/src/providers/wms/CMakeLists.txt +++ b/src/providers/wms/CMakeLists.txt @@ -4,7 +4,11 @@ SET (WMS_MOC_HDRS qgswmsprovider.h) QT4_WRAP_CPP (WMS_MOC_SRCS ${WMS_MOC_HDRS}) -INCLUDE_DIRECTORIES( . ../../core ../../core/raster ) +INCLUDE_DIRECTORIES( . + ../../core + ../../core/raster + ${GDAL_INCLUDE_DIR} +) ADD_LIBRARY(wmsprovider MODULE ${WMS_SRCS} ${WMS_MOC_SRCS})