Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
raster projector improved, also fixes #9101
  • Loading branch information
blazek committed Dec 5, 2013
1 parent fd6e37c commit 330b46b
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 35 deletions.
4 changes: 3 additions & 1 deletion src/core/raster/qgsrasterblock.cpp
Expand Up @@ -23,6 +23,8 @@
#include "qgslogger.h"
#include "qgsrasterblock.h"

const QRgb QgsRasterBlock::mNoDataColor = qRgba( 255, 255, 255, 0 );

QgsRasterBlock::QgsRasterBlock()
: mValid( true )
, mDataType( QGis::UnknownDataType )
Expand Down Expand Up @@ -464,7 +466,7 @@ bool QgsRasterBlock::setIsNoData()
return false;
}
QgsDebugMsg( "Fill image" );
mImage->fill( qRgba( 0, 0, 0, 0 ) );
mImage->fill( mNoDataColor );
return true;
}
}
Expand Down
2 changes: 2 additions & 0 deletions src/core/raster/qgsrasterblock.h
Expand Up @@ -357,6 +357,8 @@ class CORE_EXPORT QgsRasterBlock
// No data value
double mNoDataValue;

static const QRgb mNoDataColor;

// Data block for numerical data types, not used with image data types
// QByteArray does not seem to be intended for large data blocks, does it?
void * mData;
Expand Down
73 changes: 43 additions & 30 deletions src/core/raster/qgsrasterprojector.cpp
Expand Up @@ -290,6 +290,9 @@ void QgsRasterProjector::calcSrcExtent()
}
// Expand a bit to avoid possible approx coords falling out because of representation error?

// Combine with maximum source extent
mSrcExtent = mSrcExtent.intersect( &mExtent );

// If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
// align extent to src resolution to avoid jumping of reprojected pixels
// when shifting resampled grid.
Expand Down Expand Up @@ -464,19 +467,19 @@ void QgsRasterProjector::nextHelper()
mHelperTopRow++;
}

void QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
bool QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
{
if ( mApproximate )
{
approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
return approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
}
else
{
preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol, ct );
return preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol, ct );
}
}

void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
bool QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
{
#ifdef QGISDEBUG
QgsDebugMsgLevel( QString( "theDestRow = %1" ).arg( theDestRow ), 5 );
Expand All @@ -501,31 +504,31 @@ void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *
QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
#endif

if ( !mExtent.contains( QgsPoint( x, y ) ) )
{
return false;
}
// Get source row col
*theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - y ) / mSrcYRes );
*theSrcCol = ( int ) floor(( x - mSrcExtent.xMinimum() ) / mSrcXRes );
#ifdef QGISDEBUG
QgsDebugMsgLevel( QString( "mSrcExtent.yMaximum() = %1 mSrcYRes = %2" ).arg( mSrcExtent.yMaximum() ).arg( mSrcYRes ), 5 );
QgsDebugMsgLevel( QString( "mSrcExtent.yMinimum() = %1 mSrcExtent.yMaximum() = %2 mSrcYRes = %3" ).arg( mSrcExtent.yMinimum() ).arg( mSrcExtent.yMaximum() ).arg( mSrcYRes ), 5 );
QgsDebugMsgLevel( QString( "theSrcRow = %1 theSrcCol = %2" ).arg( *theSrcRow ).arg( *theSrcCol ), 5 );
#endif

// With epsg 32661 (Polar Stereographic) it was happening that *theSrcCol == mSrcCols
// For now silently correct limits to avoid crashes
// TODO: review
if ( *theSrcRow >= mSrcRows )
*theSrcRow = mSrcRows - 1;
if ( *theSrcRow < 0 )
*theSrcRow = 0;
if ( *theSrcCol >= mSrcCols )
*theSrcCol = mSrcCols - 1;
if ( *theSrcCol < 0 )
*theSrcCol = 0;
// should not happen
if ( *theSrcRow >= mSrcRows ) return false;
if ( *theSrcRow < 0 ) return false;
if ( *theSrcCol >= mSrcCols ) return false;
if ( *theSrcCol < 0 ) return false;

Q_ASSERT( *theSrcRow < mSrcRows );
Q_ASSERT( *theSrcCol < mSrcCols );
return true;
}

void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
bool QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
{
int myMatrixRow = matrixRow( theDestRow );
int myMatrixCol = matrixCol( theDestCol );
Expand Down Expand Up @@ -561,23 +564,25 @@ void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, i
double mySrcX = bx + ( tx - bx ) * yfrac;
double mySrcY = by + ( ty - by ) * yfrac;

if ( !mExtent.contains( QgsPoint( mySrcX, mySrcY ) ) )
{
return false;
}

// TODO: check again cell selection (coor is in the middle)

*theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes );
*theSrcCol = ( int ) floor(( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes );

// For now silently correct limits to avoid crashes
// TODO: review
if ( *theSrcRow >= mSrcRows )
*theSrcRow = mSrcRows - 1;
if ( *theSrcRow < 0 )
*theSrcRow = 0;
if ( *theSrcCol >= mSrcCols )
*theSrcCol = mSrcCols - 1;
if ( *theSrcCol < 0 )
*theSrcCol = 0;
Q_ASSERT( *theSrcRow < mSrcRows );
Q_ASSERT( *theSrcCol < mSrcCols );
// should not happen
if ( *theSrcRow >= mSrcRows ) return false;
if ( *theSrcRow < 0 ) return false;
if ( *theSrcCol >= mSrcCols ) return false;
if ( *theSrcCol < 0 ) return false;

return true;
}

void QgsRasterProjector::insertRows( const QgsCoordinateTransform* ct )
Expand Down Expand Up @@ -816,7 +821,11 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex
return outputBlock;
}

// No data:
// set output to no data, it should be fast
outputBlock->setIsNoData();

// No data: because isNoData()/setIsNoData() is slow with respect to simple memcpy,
// we use if only if necessary:
// 1) no data value exists (numerical) -> memcpy, not necessary isNoData()/setIsNoData()
// 2) no data value does not exist but it may contain no data (numerical no data bitmap)
// -> must use isNoData()/setIsNoData()
Expand All @@ -825,27 +834,31 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex

// To copy no data values stored in bitmaps we have to use isNoData()/setIsNoData(),
// we cannot fill output block with no data because we use memcpy for data, not setValue().
bool doNoData = inputBlock->hasNoData() && !inputBlock->hasNoDataValue();
bool doNoData = !QgsRasterBlock::typeIsNumeric( inputBlock->dataType() ) && inputBlock->hasNoData() && !inputBlock->hasNoDataValue();

const QgsCoordinateTransform* ct = 0;
if ( !mApproximate )
{
ct = 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 )
{
srcRowCol( i, j, &srcRow, &srcCol, ct );
bool inside = srcRowCol( i, j, &srcRow, &srcCol, ct );
if ( !inside ) continue; // we have everything set to no data

qgssize srcIndex = ( qgssize )srcRow * mSrcCols + 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
if ( doNoData && inputBlock->isNoData( srcRow, srcCol ) )
{
outputBlock->setIsNoData( srcRow, srcCol );
outputBlock->setIsNoData( i, j );
continue ;
}

Expand Down
11 changes: 7 additions & 4 deletions src/core/raster/qgsrasterprojector.h
Expand Up @@ -113,8 +113,11 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
void setSrcRows( int theRows ) { mSrcRows = theRows; mSrcXRes = mSrcExtent.height() / mSrcRows; }
void setSrcCols( int theCols ) { mSrcCols = theCols; mSrcYRes = mSrcExtent.width() / mSrcCols; }

/** \brief Get source row and column indexes for current source extent and resolution */
void srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct );
/** \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 );

int dstRows() const { return mDestRows; }
int dstCols() const { return mDestCols; }
Expand All @@ -130,10 +133,10 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
QgsPoint srcPoint( int theRow, int theCol );

/** \brief Get precise source row and column indexes for current source extent and resolution */
inline void preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct );
inline bool preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct );

/** \brief Get approximate source row and column indexes for current source extent and resolution */
inline void approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol );
inline bool approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol );

/** \brief Calculate matrix */
void calc();
Expand Down

0 comments on commit 330b46b

Please sign in to comment.