Skip to content

Commit

Permalink
throw away all the tricks and do second resampling, it is quite fast
Browse files Browse the repository at this point in the history 
git-svn-id: http://svn.osgeo.org/qgis/trunk@15530 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
rblazek committed Mar 17, 2011
1 parent d8a2848 commit a2020cb
Showing 1 changed file with 70 additions and 177 deletions.
247 changes: 70 additions & 177 deletions src/providers/gdal/qgsgdalprovider.cpp
View file
Expand Up @@ -38,6 +38,7 @@
#include <QFileInfo>
#include <QFile>
#include <QHash>
#include <QTime>

#include "gdalwarper.h"
#include "ogr_spatialref.h"
Expand Down Expand Up @@ -630,58 +631,26 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
double srcYRes = mGeoTransform[5]; // may be negative?
QgsDebugMsg( QString( "xRes = %1 yRes = %2 srcXRes = %3 srcYRes = %4" ).arg( xRes ).arg( yRes ).arg( srcXRes ).arg( srcYRes ) );

double center, centerRaster;

// target size in pizels
int width = right - left + 1;
int height = bottom - top + 1;

int srcWidth, srcHeight; // source size in pixels

int srcLeft = 0; // source raster x offset
int srcTop = 0; // source raster x offset
int srcBottom = ySize() - 1;
int srcRight = xSize() - 1;

// We have to extend destination grid for GDALRasterIO to keep resolution:
double topSpace = 0;
double bottomSpace = 0;
double leftSpace = 0;
double rightSpace = 0;

// It is easier to think separately about 3 cases for each axe: xRes = srcXRes, xRes > srcXRes, xRes < srcXRes, yRes = srcYRes, yRes > srcYRes, yRes < srcYRes
// Some expressions may be duplicated but it is better for keeping clear ideas
QTime time;
time.start();
// Note: original approach for xRes < srcXRes || yRes < qAbs( srcYRes ) was to avoid
// second resampling and read with GDALRasterIO to another temporary data block
// extended to fit src grid. The problem was that with src resolution much bigger
// than dst res, the target could become very large
// in theory it was going to infinity when zooming in ...

// *** CASE 1: xRes = srcXRes, yRes = srcYRes
// This is not that rare because of 'Zoom to best resolution' tool
// We just need to align the first cell
if ( xRes == srcXRes )
{
if ( mExtent.xMinimum() < myRasterExtent.xMinimum() ) // raster runs outside to left, calc offset
{
srcLeft = qRound(( myRasterExtent.xMinimum() - mExtent.xMinimum() ) / srcXRes );
}
if ( mExtent.xMaximum() > myRasterExtent.xMaximum() )
{
srcRight = qRound(( myRasterExtent.xMaximum() - mExtent.xMinimum() ) / srcXRes ) - 1 ;
}
}

if ( yRes == qAbs( srcYRes ) )
{
// GDAL states that mGeoTransform[3] is top, may it also be bottom and mGeoTransform[5] positive?
if ( mExtent.yMaximum() > myRasterExtent.yMaximum() ) // raster runs outside up, calc offset
{
//QgsDebugMsg( QString( "mExtent.yMaximum() - myRasterExtent.yMaximum() = %1 ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes = %2 srcTop = %3" ).arg( mExtent.yMaximum() - myRasterExtent.yMaximum() ).arg( ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes ).arg( srcTop ) );
srcTop = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes );
}
if ( mExtent.yMinimum() < myRasterExtent.yMinimum() )
{
srcBottom = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMinimum() ) / srcYRes ) - 1 ;
}
}

// *** CASE 2: xRes > srcXRes, yRes > srcYRes
// Note: original approach for xRes > srcXRes, yRes > srcYRes was to read directly with GDALRasterIO
// but we would face this problem:
// If the edge of the source is greater than the edge of destination:
// src: | | | | | | | | |
// dst: | | | |
Expand All @@ -695,173 +664,97 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
// case, but the first column may be rendered as without values event if its center
// is covered by src column. That could result in wrongly rendered (missing) edges
// which could be easily noticed by user
//
// Conclusion: Both are incorrect, but the b) is probably worse because it can be easily
// noticed. Missing data are worse than slightly shifted nearest src cells.
//

if ( xRes > srcXRes )
// Because of problems mentioned above we read to another temporary block and do i
// another resampling here which appeares to be quite fast

// Get necessary src extent aligned to src resolution
if ( mExtent.xMinimum() < myRasterExtent.xMinimum() )
{
if ( mExtent.xMinimum() < myRasterExtent.xMinimum() )
{
srcLeft = qRound(( myRasterExtent.xMinimum() - mExtent.xMinimum() ) / srcXRes );
}
if ( mExtent.xMaximum() > myRasterExtent.xMaximum() )
{
srcRight = qRound(( myRasterExtent.xMaximum() - mExtent.xMinimum() ) / srcXRes ) - 1 ;
}
srcLeft = static_cast<int>( floor(( myRasterExtent.xMinimum() - mExtent.xMinimum() ) / srcXRes ) );
}
if ( yRes > qAbs( srcYRes ) )
if ( mExtent.xMaximum() > myRasterExtent.xMaximum() )
{
if ( mExtent.yMaximum() > myRasterExtent.yMaximum() )
{
srcTop = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes );
}
if ( mExtent.yMinimum() < myRasterExtent.yMinimum() )
{
srcBottom = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMinimum() ) / srcYRes ) - 1 ;
}
srcRight = static_cast<int>( floor(( myRasterExtent.xMaximum() - mExtent.xMinimum() ) / srcXRes ) );
}

// *** CASE 3: xRes < srcXRes, yRes < srcYRes
// IMHO, we cannot align target grid to raster grid using target grid edges
// and finding the nearest raster grid because it could lead to cell center
// getting outside the right cell when doing resampling, example
// src: | | | |
// dst: | | | |
// Raster width is 30m and it has 3 columns and we want to read xrange 5.1-30
// to 3 columns, the nearest edge for beginning in raster grid is 10.0
// reading cols 1-2, we get raster[1] value in target[0], but the center of
// target[0] is 5.1 + ((30-5.1)/3)/2 = 9.25 so it falls to raster[0]. Is it right?
// => We are looking for such alignment with which the center of first/last cell
// falls to the right raster cell
int xAddPixels = 0;
int leftAddPixels = 0;
if ( xRes < srcXRes )
// GDAL states that mGeoTransform[3] is top, may it also be bottom and mGeoTransform[5] positive?
if ( mExtent.yMaximum() > myRasterExtent.yMaximum() )
{
center = myRasterExtent.xMinimum() + xRes / 2;
centerRaster = ( center - mExtent.xMinimum() ) / srcXRes;
srcLeft = static_cast<int>( floor( centerRaster ) );

// Space which has to be added to destination to fit better source
// Warning - TODO: if xRes is much smaller than srcXRes, adding space and pixels to dst may
// result in very large dst grids!!! E.g.
// src: | | |
// dst: | | |
// Solution could be vector rendering of rectangles.

leftSpace = myRasterExtent.xMinimum() - ( mExtent.xMinimum() + srcLeft * srcXRes );
leftSpace = leftSpace > 0 ? leftSpace : 0; // makes only sense for positive
center = myRasterExtent.xMaximum() - xRes / 2;
centerRaster = ( center - mExtent.xMinimum() ) / srcXRes;
srcRight = static_cast<int>( floor( centerRaster ) );

rightSpace = ( mExtent.xMinimum() + ( srcRight + 1 ) * srcXRes ) - myRasterExtent.xMaximum();
rightSpace = rightSpace > 0 ? rightSpace : 0;
QgsDebugMsg( QString( "center = %1 centerRaster = %2 srcRight = %3 rightSpace = %4" ).arg( center ).arg( centerRaster ).arg( srcRight ).arg( rightSpace ) );

double xAdd = leftSpace + rightSpace;
xAddPixels = qRound( xAdd / xRes );
leftAddPixels = qRound( leftSpace / xRes );
srcTop = static_cast<int>( floor( -1. * ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes ) );
}

int yAddPixels = 0;
int topAddPixels = 0;
if ( yRes < qAbs( srcYRes ) )
if ( mExtent.yMinimum() < myRasterExtent.yMinimum() )
{
center = myRasterExtent.yMaximum() - yRes / 2;
centerRaster = -1. * ( mExtent.yMaximum() - center ) / srcYRes;
srcTop = static_cast<int>( floor( centerRaster ) );
topSpace = ( mExtent.yMaximum() + srcTop * srcYRes ) - myRasterExtent.yMaximum();
topSpace = topSpace > 0 ? topSpace : 0;
center = myRasterExtent.yMinimum() + yRes / 2;
centerRaster = -1. * ( mExtent.yMaximum() - center ) / srcYRes;

srcBottom = static_cast<int>( floor( centerRaster ) );
QgsDebugMsg( QString( "myRasterExtent.yMinimum() = %1 myRasterExtent.yMaximum() = %2" ).arg( myRasterExtent.yMinimum() ).arg( myRasterExtent.yMaximum() ) );
bottomSpace = myRasterExtent.yMinimum() - ( mExtent.yMaximum() + ( srcBottom + 1 ) * srcYRes );
bottomSpace = bottomSpace > 0 ? bottomSpace : 0;
QgsDebugMsg( QString( "center = %1 centerRaster = %2 srcBottom = %3 bottomSpace = %4" ).arg( center ).arg( centerRaster ).arg( srcBottom ).arg( bottomSpace ) );

double yAdd = topSpace + bottomSpace;
yAddPixels = qRound( yAdd / yRes );
topAddPixels = qRound( topSpace / yRes );
srcBottom = static_cast<int>( floor( -1. * ( mExtent.yMaximum() - myRasterExtent.yMinimum() ) / srcYRes ) );
}

srcWidth = srcRight - srcLeft + 1;
srcHeight = srcBottom - srcTop + 1;

QgsDebugMsg( QString( "srcTop = %1 srcBottom = %2 srcLeft = %3 srcRight = %4" ).arg( srcTop ).arg( srcBottom ).arg( srcLeft ).arg( srcRight ) );

QgsDebugMsg( QString( "topSpace = %1 bottomSpace = %2 leftSpace = %3 rightSpace = %4" ).arg( topSpace ).arg( bottomSpace ).arg( leftSpace ).arg( rightSpace ) );

int srcWidth = srcRight - srcLeft + 1;
int srcHeight = srcBottom - srcTop + 1;

QgsDebugMsg( QString( "width = %1 height = %2 srcWidth = %3 srcHeight = %4" ).arg( width ).arg( height ).arg( srcWidth ).arg( srcHeight ) );

int tmpWidth = srcWidth;
int tmpHeight = srcHeight;

QgsDebugMsg( QString( "xAddPixels = %1 yAddPixels = %2 leftAddPixels = %3 topAddPixels = %4" ).arg( xAddPixels ).arg( yAddPixels ).arg( leftAddPixels ).arg( topAddPixels ) );

int totalWidth = width + xAddPixels;
int totalHeight = height + yAddPixels;

QgsDebugMsg( QString( "totalWidth = %1 totalHeight = %2" ).arg( totalWidth ).arg( totalHeight ) );
if ( xRes > srcXRes ) tmpWidth = width;
if ( yRes > srcYRes ) tmpHeight = height;

// Allocate temporary block
char *tmpBlock = ( char * )malloc( dataSize * tmpWidth * tmpHeight );

GDALRasterBandH gdalBand = GDALGetRasterBand( mGdalDataset, theBandNo );
GDALDataType type = ( GDALDataType )mGdalDataType[theBandNo-1];
CPLErrorReset();
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
srcLeft, srcTop, srcWidth, srcHeight,
( void * )tmpBlock,
tmpWidth, tmpHeight, type,
0, 0 );

if ( xAddPixels == 0 && yAddPixels == 0 )
if ( err != CPLE_None )
{
// Calc beginnig of data if raster does not start at top
block = ( char * ) theBlock;
if ( top != 0 )
{
block += dataSize * thePixelWidth * top;
}

// Cal nLineSpace if raster does not cover whole extent
int nLineSpace = dataSize * thePixelWidth;
if ( left != 0 )
{
block += dataSize * left;
}
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
srcLeft, srcTop, srcWidth, srcHeight,
( void * )block,
width, height, type,
0, nLineSpace );
QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
QgsDebugMsg( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
free( tmpBlock );
return;
}
else
{
// Allocate temporary block
void *tmpBlock = malloc( dataSize * totalWidth * totalHeight );

CPLErrorReset();
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
srcLeft, srcTop, srcWidth, srcHeight,
( void * )tmpBlock,
totalWidth, totalHeight, type,
0, 0 );
QgsDebugMsg( QString( "GDALRasterIO time (ms): %1" ).arg( time.elapsed() ) );
time.start();

if ( err != CPLE_None )
{
QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
QgsDebugMsg( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
free( tmpBlock );
return;
}
double tmpXMin = mExtent.xMinimum() + srcLeft * srcXRes;
double tmpYMax = mExtent.yMaximum() + srcTop * srcYRes;
double tmpXRes = srcWidth * srcXRes / tmpWidth;
double tmpYRes = srcHeight * srcYRes / tmpHeight; // negative

QgsDebugMsg( QString( "tmpXMin = %1 tmpYMax = %2 tmpWidth = %3 tmpHeight = %4" ).arg( tmpXMin ).arg( tmpYMax ).arg( tmpWidth ).arg( tmpHeight ) );

for ( int i = 0; i < height; i++ )
for ( int row = 0; row < height; row++ )
{
double y = myRasterExtent.yMaximum() - ( row + 0.5 ) * yRes;
int tmpRow = static_cast<int>( floor( -1. * ( tmpYMax - y ) / tmpYRes ) );

char *srcRowBlock = tmpBlock + dataSize * tmpRow * tmpWidth;
char *dstRowBlock = ( char * )theBlock + dataSize * ( top + row ) * thePixelWidth;
for ( int col = 0; col < width; col++ )
{
int r = i + topAddPixels;
char *src = ( char * )tmpBlock + dataSize * r * totalWidth + dataSize * leftAddPixels;
char *dst = ( char * )theBlock + dataSize * ( top + i ) * thePixelWidth + dataSize * ( left );
memcpy( dst, src, dataSize*width );
// cell center
double x = myRasterExtent.xMinimum() + ( col + 0.5 ) * xRes;
// floor() is quite slow! Use just cast to int.
int tmpCol = static_cast<int>(( x - tmpXMin ) / tmpXRes ) ;
//QgsDebugMsg( QString( "row = %1 col = %2 tmpRow = %3 tmpCol = %4" ).arg(row).arg(col).arg(tmpRow).arg(tmpCol) );

char *src = srcRowBlock + dataSize * tmpCol;
char *dst = dstRowBlock + dataSize * ( left + col );
memcpy( dst, src, dataSize );
}

free( tmpBlock );
}

free( tmpBlock );
QgsDebugMsg( QString( "resample time (ms): %1" ).arg( time.elapsed() ) );

return;
}

Expand Down

0 comments on commit a2020cb

Please sign in to comment.