From fdeac8198c91c2a15385c124b8298c14677ca431 Mon Sep 17 00:00:00 2001 From: Martin Dobias Date: Thu, 1 Sep 2016 18:25:01 +0800 Subject: [PATCH] Fix incorrect raster reprojection after drawing raster previews Moved all temporary projector members to a private class, so even recursive block() calls will not affect each other (there is no new performance penalty as block() call always recomputes the temporary control point matrix anyway) --- python/core/raster/qgsrasterprojector.sip | 3 - src/core/raster/qgsrasterprojector.cpp | 151 ++++++++++------------ src/core/raster/qgsrasterprojector.h | 74 ++++++----- 3 files changed, 108 insertions(+), 120 deletions(-) diff --git a/python/core/raster/qgsrasterprojector.sip b/python/core/raster/qgsrasterprojector.sip index 813eadfbbf52..bbc23587be34 100644 --- a/python/core/raster/qgsrasterprojector.sip +++ b/python/core/raster/qgsrasterprojector.sip @@ -43,9 +43,6 @@ class QgsRasterProjector : QgsRasterInterface /** \brief Get destination CRS */ QgsCoordinateReferenceSystem destinationCrs() const; - /** \brief set maximum source resolution */ - void setMaxSrcRes( double theMaxSrcXRes, double theMaxSrcYRes ); - Precision precision() const; void setPrecision( Precision precision ); // Translated precision mode, for use in ComboBox etc. diff --git a/src/core/raster/qgsrasterprojector.cpp b/src/core/raster/qgsrasterprojector.cpp index 626333b6f8d7..37315926ad6b 100644 --- a/src/core/raster/qgsrasterprojector.cpp +++ b/src/core/raster/qgsrasterprojector.cpp @@ -24,31 +24,11 @@ #include "qgscsexception.h" - QgsRasterProjector::QgsRasterProjector() : QgsRasterInterface( nullptr ) , mSrcDatumTransform( -1 ) , mDestDatumTransform( -1 ) , mPrecision( Approximate ) - , mApproximate( false ) - , mDestRows( 0 ) - , mDestCols( 0 ) - , mDestXRes( 0.0 ) - , mDestYRes( 0.0 ) - , mSrcRows( 0 ) - , mSrcCols( 0 ) - , mSrcXRes( 0.0 ) - , mSrcYRes( 0.0 ) - , mDestRowsPerMatrixRow( 0.0 ) - , mDestColsPerMatrixCol( 0.0 ) - , pHelperTop( nullptr ) - , pHelperBottom( nullptr ) - , mHelperTopRow( 0 ) - , mCPCols( 0 ) - , mCPRows( 0 ) - , mSqrTolerance( 0.0 ) - , mMaxSrcXRes( 0 ) - , mMaxSrcYRes( 0 ) { QgsDebugMsgLevel( "Entered", 4 ); } @@ -62,17 +42,12 @@ QgsRasterProjector* QgsRasterProjector::clone() const projector->mDestCRS = mDestCRS; projector->mSrcDatumTransform = mSrcDatumTransform; projector->mDestDatumTransform = mDestDatumTransform; - projector->mMaxSrcXRes = mMaxSrcXRes; - projector->mMaxSrcYRes = mMaxSrcYRes; - projector->mExtent = mExtent; projector->mPrecision = mPrecision; return projector; } QgsRasterProjector::~QgsRasterProjector() { - delete[] pHelperTop; - delete[] pHelperBottom; } int QgsRasterProjector::bandCount() const @@ -97,22 +72,36 @@ void QgsRasterProjector::setCrs( const QgsCoordinateReferenceSystem & theSrcCRS, mDestDatumTransform = destDatumTransform; } -void QgsRasterProjector::calc() + +ProjectorData::ProjectorData( const QgsRectangle& extent, int width, int height, QgsRasterInterface* input, const QgsCoordinateTransform& inverseCt, QgsRasterProjector::Precision precision ) + : mApproximate( false ) + , mInverseCt( new QgsCoordinateTransform( inverseCt ) ) + , mDestExtent( extent ) + , mDestRows( height ) + , mDestCols( width ) + , mDestXRes( 0.0 ) + , mDestYRes( 0.0 ) + , mSrcRows( 0 ) + , mSrcCols( 0 ) + , mSrcXRes( 0.0 ) + , mSrcYRes( 0.0 ) + , mDestRowsPerMatrixRow( 0.0 ) + , mDestColsPerMatrixCol( 0.0 ) + , pHelperTop( nullptr ) + , pHelperBottom( nullptr ) + , mHelperTopRow( 0 ) + , mCPCols( 0 ) + , mCPRows( 0 ) + , mSqrTolerance( 0.0 ) + , mMaxSrcXRes( 0 ) + , mMaxSrcYRes( 0 ) { QgsDebugMsgLevel( "Entered", 4 ); - mCPMatrix.clear(); - mCPLegalMatrix.clear(); - delete[] pHelperTop; - pHelperTop = nullptr; - delete[] pHelperBottom; - pHelperBottom = nullptr; // Get max source resolution and extent if possible - mMaxSrcXRes = 0; - mMaxSrcYRes = 0; - if ( mInput ) + if ( input ) { - QgsRasterDataProvider *provider = dynamic_cast( mInput->sourceInput() ); + QgsRasterDataProvider *provider = dynamic_cast( input->sourceInput() ); if ( provider ) { if ( provider->capabilities() & QgsRasterDataProvider::Size ) @@ -138,9 +127,7 @@ void QgsRasterProjector::calc() double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes; mSqrTolerance = myDestRes * myDestRes; - QgsCoordinateTransform inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform ); - - if ( mPrecision == Approximate ) + if ( precision == QgsRasterProjector::Approximate ) { mApproximate = true; } @@ -219,7 +206,15 @@ void QgsRasterProjector::calc() mSrcXRes = mSrcExtent.width() / mSrcCols; } -void QgsRasterProjector::calcSrcExtent() +ProjectorData::~ProjectorData() +{ + delete[] pHelperTop; + delete[] pHelperBottom; + delete mInverseCt; +} + + +void ProjectorData::calcSrcExtent() { /* Run around the mCPMatrix and find source extent */ // Attention, source limits are not necessarily on destination edges, e.g. @@ -283,7 +278,7 @@ void QgsRasterProjector::calcSrcExtent() QgsDebugMsgLevel( "mSrcExtent = " + mSrcExtent.toString(), 4 ); } -QString QgsRasterProjector::cpToString() +QString ProjectorData::cpToString() { QString myString; for ( int i = 0; i < mCPRows; i++ ) @@ -308,7 +303,7 @@ QString QgsRasterProjector::cpToString() return myString; } -void QgsRasterProjector::calcSrcRowsCols() +void ProjectorData::calcSrcRowsCols() { // Wee need to calculate minimum cell size in the source // TODO: Think it over better, what is the right source resolution? @@ -347,11 +342,10 @@ void QgsRasterProjector::calcSrcRowsCols() else { // take highest from corners, points in in the middle of corners and center (3 x 3 ) - QgsCoordinateTransform inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform ); //double QgsRectangle srcExtent; int srcXSize, srcYSize; - if ( extentSize( inverseCt, mDestExtent, mDestCols, mDestRows, srcExtent, srcXSize, srcYSize ) ) + if ( QgsRasterProjector::extentSize( *mInverseCt, mDestExtent, mDestCols, mDestRows, srcExtent, srcXSize, srcYSize ) ) { double srcXRes = srcExtent.width() / srcXSize; double srcYRes = srcExtent.height() / srcYSize; @@ -383,29 +377,22 @@ void QgsRasterProjector::calcSrcRowsCols() } -inline void QgsRasterProjector::destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY ) +inline void ProjectorData::destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY ) { *theX = mDestExtent.xMinimum() + theCol * mDestExtent.width() / ( mCPCols - 1 ); *theY = mDestExtent.yMaximum() - theRow * mDestExtent.height() / ( mCPRows - 1 ); } -inline int QgsRasterProjector::matrixRow( int theDestRow ) +inline int ProjectorData::matrixRow( int theDestRow ) { return static_cast< int >( floor(( theDestRow + 0.5 ) / mDestRowsPerMatrixRow ) ); } -inline int QgsRasterProjector::matrixCol( int theDestCol ) +inline int ProjectorData::matrixCol( int theDestCol ) { return static_cast< int >( floor(( theDestCol + 0.5 ) / mDestColsPerMatrixCol ) ); } -QgsPoint QgsRasterProjector::srcPoint( int theDestRow, int theCol ) -{ - Q_UNUSED( theDestRow ); - Q_UNUSED( theCol ); - return QgsPoint(); -} - -void QgsRasterProjector::calcHelper( int theMatrixRow, QgsPoint *thePoints ) +void ProjectorData::calcHelper( int theMatrixRow, QgsPoint *thePoints ) { // TODO?: should we also precalc dest cell center coordinates for x and y? for ( int myDestCol = 0; myDestCol < mDestCols; myDestCol++ ) @@ -430,7 +417,8 @@ void QgsRasterProjector::calcHelper( int theMatrixRow, QgsPoint *thePoints ) thePoints[myDestCol].setY( t ); } } -void QgsRasterProjector::nextHelper() + +void ProjectorData::nextHelper() { // We just switch pHelperTop and pHelperBottom, memory is not lost QgsPoint *tmp; @@ -441,7 +429,7 @@ void QgsRasterProjector::nextHelper() mHelperTopRow++; } -bool QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform& ct ) +bool ProjectorData::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ) { if ( mApproximate ) { @@ -449,11 +437,11 @@ bool QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcR } else { - return preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol, ct ); + return preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol ); } } -bool QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform& ct ) +bool ProjectorData::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ) { #ifdef QGISDEBUG QgsDebugMsgLevel( QString( "theDestRow = %1" ).arg( theDestRow ), 5 ); @@ -469,9 +457,9 @@ bool QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int * QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 ); #endif - if ( ct.isValid() ) + if ( mInverseCt->isValid() ) { - ct.transformInPlace( x, y, z ); + mInverseCt->transformInPlace( x, y, z ); } #ifdef QGISDEBUG @@ -502,7 +490,7 @@ bool QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int * return true; } -bool QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ) +bool ProjectorData::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ) { int myMatrixRow = matrixRow( theDestRow ); int myMatrixCol = matrixCol( theDestCol ); @@ -559,7 +547,7 @@ bool QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, i return true; } -void QgsRasterProjector::insertRows( const QgsCoordinateTransform& ct ) +void ProjectorData::insertRows( const QgsCoordinateTransform& ct ) { for ( int r = 0; r < mCPRows - 1; r++ ) { @@ -583,7 +571,7 @@ void QgsRasterProjector::insertRows( const QgsCoordinateTransform& ct ) } } -void QgsRasterProjector::insertCols( const QgsCoordinateTransform& ct ) +void ProjectorData::insertCols( const QgsCoordinateTransform& ct ) { for ( int r = 0; r < mCPRows; r++ ) { @@ -601,7 +589,7 @@ void QgsRasterProjector::insertCols( const QgsCoordinateTransform& ct ) } -void QgsRasterProjector::calcCP( int theRow, int theCol, const QgsCoordinateTransform& ct ) +void ProjectorData::calcCP( int theRow, int theCol, const QgsCoordinateTransform& ct ) { double myDestX, myDestY; destPointOnCPMatrix( theRow, theCol, &myDestX, &myDestY ); @@ -626,7 +614,7 @@ void QgsRasterProjector::calcCP( int theRow, int theCol, const QgsCoordinateTran } } -bool QgsRasterProjector::calcRow( int theRow, const QgsCoordinateTransform& ct ) +bool ProjectorData::calcRow( int theRow, const QgsCoordinateTransform& ct ) { QgsDebugMsgLevel( QString( "theRow = %1" ).arg( theRow ), 3 ); for ( int i = 0; i < mCPCols; i++ ) @@ -637,7 +625,7 @@ bool QgsRasterProjector::calcRow( int theRow, const QgsCoordinateTransform& ct ) return true; } -bool QgsRasterProjector::calcCol( int theCol, const QgsCoordinateTransform& ct ) +bool ProjectorData::calcCol( int theCol, const QgsCoordinateTransform& ct ) { QgsDebugMsgLevel( QString( "theCol = %1" ).arg( theCol ), 3 ); for ( int i = 0; i < mCPRows; i++ ) @@ -648,7 +636,7 @@ bool QgsRasterProjector::calcCol( int theCol, const QgsCoordinateTransform& ct ) return true; } -bool QgsRasterProjector::checkCols( const QgsCoordinateTransform& ct ) +bool ProjectorData::checkCols( const QgsCoordinateTransform& ct ) { if ( !ct.isValid() ) { @@ -693,7 +681,7 @@ bool QgsRasterProjector::checkCols( const QgsCoordinateTransform& ct ) return true; } -bool QgsRasterProjector::checkRows( const QgsCoordinateTransform& ct ) +bool ProjectorData::checkRows( const QgsCoordinateTransform& ct ) { if ( !ct.isValid() ) { @@ -766,22 +754,21 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex return mInput->block( bandNo, extent, width, height, feedback ); } - mDestExtent = extent; - mDestRows = height; - mDestCols = width; - calc(); + QgsCoordinateTransform inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform ); + + ProjectorData pd( extent, width, height, mInput, inverseCt, mPrecision ); - QgsDebugMsgLevel( QString( "srcExtent:\n%1" ).arg( mSrcExtent.toString() ), 4 ); - QgsDebugMsgLevel( QString( "srcCols = %1 srcRows = %2" ).arg( mSrcCols ).arg( mSrcRows ), 4 ); + QgsDebugMsgLevel( QString( "srcExtent:\n%1" ).arg( pd.srcExtent().toString() ), 4 ); + QgsDebugMsgLevel( QString( "srcCols = %1 srcRows = %2" ).arg( pd.srcCols() ).arg( pd.srcRows() ), 4 ); // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers - if ( mSrcRows <= 0 || mSrcCols <= 0 ) + if ( pd.srcRows() <= 0 || pd.srcCols() <= 0 ) { QgsDebugMsgLevel( "Zero srcRows or srcCols", 4 ); return new QgsRasterBlock(); } - QgsRasterBlock *inputBlock = mInput->block( bandNo, mSrcExtent, mSrcCols, mSrcRows, feedback ); + QgsRasterBlock *inputBlock = mInput->block( bandNo, pd.srcExtent(), pd.srcCols(), pd.srcRows(), feedback ); if ( !inputBlock || inputBlock->isEmpty() ) { QgsDebugMsg( "No raster data!" ); @@ -822,12 +809,6 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex // we cannot fill output block with no data because we use memcpy for data, not setValue(). bool doNoData = !QgsRasterBlock::typeIsNumeric( inputBlock->dataType() ) && inputBlock->hasNoData() && !inputBlock->hasNoDataValue(); - QgsCoordinateTransform inverseCt; - if ( !mApproximate ) - { - inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform ); - } - outputBlock->setIsNoData(); int srcRow, srcCol; @@ -835,10 +816,10 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex { for ( int j = 0; j < width; ++j ) { - bool inside = srcRowCol( i, j, &srcRow, &srcCol, inverseCt ); + bool inside = pd.srcRowCol( i, j, &srcRow, &srcCol ); if ( !inside ) continue; // we have everything set to no data - qgssize srcIndex = static_cast< qgssize >( srcRow ) * mSrcCols + srcCol; + qgssize srcIndex = static_cast< qgssize >( srcRow ) * pd.srcCols() + srcCol; QgsDebugMsgLevel( QString( "row = %1 col = %2 srcRow = %3 srcCol = %4" ).arg( i ).arg( j ).arg( srcRow ).arg( srcCol ), 5 ); // isNoData() may be slow so we check doNoData first diff --git a/src/core/raster/qgsrasterprojector.h b/src/core/raster/qgsrasterprojector.h index 9682d758ce49..b7832a68918d 100644 --- a/src/core/raster/qgsrasterprojector.h +++ b/src/core/raster/qgsrasterprojector.h @@ -74,13 +74,6 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface /** \brief Get destination CRS */ QgsCoordinateReferenceSystem destinationCrs() const { return mDestCRS; } - /** \brief set maximum source resolution */ - void setMaxSrcRes( double theMaxSrcXRes, double theMaxSrcYRes ) - { - mMaxSrcXRes = theMaxSrcXRes; - mMaxSrcYRes = theMaxSrcYRes; - } - Precision precision() const { return mPrecision; } void setPrecision( Precision precision ) { mPrecision = precision; } // Translated precision mode, for use in ComboBox etc. @@ -99,12 +92,48 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface private: + /** Source CRS */ + QgsCoordinateReferenceSystem mSrcCRS; + + /** Destination CRS */ + QgsCoordinateReferenceSystem mDestCRS; + + /** Source datum transformation id (or -1 if none) */ + int mSrcDatumTransform; + + /** Destination datum transformation id (or -1 if none) */ + int mDestDatumTransform; + + /** Requested precision */ + Precision mPrecision; + +}; + +/// @cond PRIVATE + +/** + * Internal class for reprojection of rasters - either exact or approximate. + * QgsRasterProjector creates it and then keeps calling srcRowCol() to get source pixel position + * for every destination pixel position. + */ +class ProjectorData +{ + public: + /** Initialize reprojector and calculate matrix */ + ProjectorData( const QgsRectangle &extent, int width, int height, QgsRasterInterface *input, const QgsCoordinateTransform &inverseCt, QgsRasterProjector::Precision precision ); + ~ProjectorData(); + /** \brief Get source row and column indexes for current source extent and resolution If source pixel is outside source extent theSrcRow and theSrcCol are left unchanged. @return true if inside source */ - bool srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform& ct ); + bool srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ); + + QgsRectangle srcExtent() const { return mSrcExtent; } + int srcRows() const { return mSrcRows; } + int srcCols() const { return mSrcCols; } + private: /** \brief get destination point for _current_ destination position */ void destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY ); @@ -112,18 +141,12 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface int matrixRow( int theDestRow ); int matrixCol( int theDestCol ); - /** \brief get destination point for _current_ matrix position */ - QgsPoint srcPoint( int theRow, int theCol ); - /** \brief Get precise source row and column indexes for current source extent and resolution */ - inline bool preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform& ct ); + inline bool preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ); /** \brief Get approximate source row and column indexes for current source extent and resolution */ inline bool approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol ); - /** \brief Calculate matrix */ - void calc(); - /** \brief insert rows to matrix */ void insertRows( const QgsCoordinateTransform& ct ); @@ -162,27 +185,12 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface /** Get mCPMatrix as string */ QString cpToString(); - /** Source CRS */ - QgsCoordinateReferenceSystem mSrcCRS; - - /** Destination CRS */ - QgsCoordinateReferenceSystem mDestCRS; - - /** Source datum transformation id (or -1 if none) */ - int mSrcDatumTransform; - - /** Destination datum transformation id (or -1 if none) */ - int mDestDatumTransform; - - /** Requested precision */ - Precision mPrecision; - /** Use approximation (requested precision is Approximate and it is possible to calculate * an approximation matrix with a sufficient precision) */ bool mApproximate; - - + /** Transformation from destination CRS to source CRS */ + QgsCoordinateTransform* mInverseCt; /** Destination extent */ QgsRectangle mDestExtent; @@ -255,5 +263,7 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface }; +/// @endcond + #endif