Skip to content
Permalink
Browse files
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)
  • Loading branch information
@wonder-sk
wonder-sk committed Sep 1, 2016
1 parent 6bb291f commit fdeac81
Show file tree
Hide file tree
Showing 3 changed files with 108 additions and 120 deletions.
3 python/core/raster/qgsrasterprojector.sip
View file
@@ -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.
151 src/core/raster/qgsrasterprojector.cpp
View file
@@ -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<QgsRasterDataProvider*>( mInput->sourceInput() );
QgsRasterDataProvider *provider = dynamic_cast<QgsRasterDataProvider*>( 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,19 +429,19 @@ 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 )
{
return approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
}
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,23 +809,17 @@ 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;
for ( int i = 0; i < height; ++i )
{
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

0 comments on commit fdeac81

Please sign in to comment.