Skip to content

Commit c4ccb8d

Browse files
author
rblazek
committedMar 16, 2011
better resampling alignment
git-svn-id: http://svn.osgeo.org/qgis/trunk@15514 c8812cc2-4d05-0410-92ff-de0c093fc19c
·
ltr-3_40Before-merge-new_vector_api
1 parent a9698d8 commit c4ccb8d

File tree

1 file changed

+210
-123
lines changed

1 file changed

+210
-123
lines changed
 

‎src/providers/gdal/qgsgdalprovider.cpp

Lines changed: 210 additions & 123 deletions
Original file line numberDiff line numberDiff line change
@@ -348,13 +348,13 @@ bool QgsGdalProvider::crsFromWkt( const char *wkt )
348348
{
349349
void *hCRS = OSRNewSpatialReference( NULL );
350350

351-
if ( OSRImportFromWkt( hCRS, (char **) &wkt ) == OGRERR_NONE )
351+
if ( OSRImportFromWkt( hCRS, ( char ** ) &wkt ) == OGRERR_NONE )
352352
{
353353
if ( OSRAutoIdentifyEPSG( hCRS ) == OGRERR_NONE )
354354
{
355355
QString authid = QString( "%1:%2" )
356-
.arg( OSRGetAuthorityName( hCRS, NULL ) )
357-
.arg( OSRGetAuthorityCode( hCRS, NULL ) );
356+
.arg( OSRGetAuthorityName( hCRS, NULL ) )
357+
.arg( OSRGetAuthorityCode( hCRS, NULL ) );
358358
QgsDebugMsg( "authid recognized as " + authid );
359359
mCrs.createFromOgcWmsCrs( authid );
360360
}
@@ -570,14 +570,25 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
570570
QgsDebugMsg( QString( "transform : %1" ).arg( mGeoTransform[i] ) );
571571
}
572572

573-
// TODO: fill block with no data values
573+
int dataSize = dataTypeSize( theBandNo ) / 8;
574+
575+
// fill with null values
576+
QByteArray ba = noValueBytes( theBandNo );
577+
char *nodata = ba.data();
578+
char *block = ( char * ) theBlock;
579+
for ( int i = 0; i < thePixelWidth * thePixelHeight; i++ )
580+
{
581+
memcpy( block, nodata, dataSize );
582+
block += dataSize;
583+
}
574584

575585
QgsRectangle myRasterExtent = theExtent.intersect( &mExtent );
576586
if ( myRasterExtent.isEmpty() )
577587
{
578588
QgsDebugMsg( "draw request outside view extent." );
579589
return;
580590
}
591+
QgsDebugMsg( "mExtent: " + mExtent.toString() );
581592
QgsDebugMsg( "myRasterExtent: " + myRasterExtent.toString() );
582593

583594
double xRes = theExtent.width() / thePixelWidth;
@@ -614,167 +625,243 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
614625

615626
// Calculate rows/cols limits in raster grid space
616627

617-
// IMHO, we cannot align target grid to raster grid using target grid edges
618-
// and finding the nearest raster grid because it could lead to cell center
619-
// getting outside the right cell when doing resampling, example
620-
// Raster width is 30m and it has 3 columns and we want to read xrange 5.1-30
621-
// to 3 columns, the nearest edge for beginning in raster grid is 10.0
622-
// reading cols 1-2, we get raster[1] value in target[0], but the center of
623-
// target[0] is 5.1 + ((30-5.1)/3)/2 = 9.25 so it falls to raster[0]. Is it right?
624-
// => We are looking for such alignment with which the center of first/last cell
625-
// alls to the right raster cell
628+
// Set readable names
629+
double srcXRes = mGeoTransform[1];
630+
double srcYRes = mGeoTransform[5]; // may be negative?
631+
QgsDebugMsg( QString( "xRes = %1 yRes = %2 srcXRes = %3 srcYRes = %4" ).arg( xRes ).arg( yRes ).arg( srcXRes ).arg( srcYRes ) );
626632

627633
double center, centerRaster;
628-
int rasterTop, rasterBottom, rasterLeft, rasterRight;
634+
635+
// target size in pizels
636+
int width = right - left + 1;
637+
int height = bottom - top + 1;
638+
639+
int srcWidth, srcHeight; // source size in pixels
640+
641+
int srcLeft = 0; // source raster x offset
642+
int srcTop = 0; // source raster x offset
643+
int srcBottom = ySize() - 1;
644+
int srcRight = xSize() - 1;
629645

630646
// We have to extend destination grid for GDALRasterIO to keep resolution:
631-
double topSpace, bottomSpace, leftSpace, rightSpace;
632-
633-
// top
634-
center = myRasterExtent.yMaximum() - yRes / 2;
635-
// center in raster space
636-
// GDAL states that mGeoTransform[3] is top, but probably it can be also bottom
637-
// if mGeoTransform[5] is negative ??? No, mGeoTransform[5] is negative normaly
638-
// - vice versa?
639-
//if ( mGeoTransform[5] > 0 )
640-
if ( mGeoTransform[5] < 0 )
647+
double topSpace = 0;
648+
double bottomSpace = 0;
649+
double leftSpace = 0;
650+
double rightSpace = 0;
651+
652+
// It is easier to think separately about 3 cases for each axe: xRes = srcXRes, xRes > srcXRes, xRes < srcXRes, yRes = srcYRes, yRes > srcYRes, yRes < srcYRes
653+
// Some expressions may be duplicated but it is better for keeping clear ideas
654+
655+
// *** CASE 1: xRes = srcXRes, yRes = srcYRes
656+
// This is not that rare because of 'Zoom to best resolution' tool
657+
// We just need to align the first cell
658+
if ( xRes == srcXRes )
641659
{
642-
centerRaster = -1. * ( mGeoTransform[3] - center ) / mGeoTransform[5];
660+
if ( mExtent.xMinimum() < myRasterExtent.xMinimum() ) // raster runs outside to left, calc offset
661+
{
662+
srcLeft = qRound(( myRasterExtent.xMinimum() - mExtent.xMinimum() ) / srcXRes );
663+
}
664+
if ( mExtent.xMaximum() > myRasterExtent.xMaximum() )
665+
{
666+
srcRight = qRound(( myRasterExtent.xMaximum() - mExtent.xMinimum() ) / srcXRes ) - 1 ;
667+
}
643668
}
644-
else
669+
670+
if ( yRes == qAbs( srcYRes ) )
645671
{
646-
centerRaster = ( center - mGeoTransform[3] ) / mGeoTransform[5];
672+
// GDAL states that mGeoTransform[3] is top, may it also be bottom and mGeoTransform[5] positive?
673+
if ( mExtent.yMaximum() > myRasterExtent.yMaximum() ) // raster runs outside up, calc offset
674+
{
675+
//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 ) );
676+
srcTop = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes );
677+
}
678+
if ( mExtent.yMinimum() < myRasterExtent.yMinimum() )
679+
{
680+
srcBottom = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMinimum() ) / srcYRes ) - 1 ;
681+
}
647682
}
648-
rasterTop = static_cast<int>( floor( centerRaster ) );
649-
topSpace = ( mGeoTransform[3] + rasterTop * mGeoTransform[5] ) - myRasterExtent.yMaximum();
650683

651-
// bottom
652-
center = myRasterExtent.yMinimum() + yRes / 2;
653-
//if ( mGeoTransform[5] > 0 )
654-
if ( mGeoTransform[5] < 0 )
684+
// *** CASE 2: xRes > srcXRes, yRes > srcYRes
685+
// If the edge of the source is greater than the edge of destination:
686+
// src: | | | | | | | | |
687+
// dst: | | | |
688+
// We have 2 options for resampling:
689+
// a) 'Stretch' the src and align the start edge of src to the start edge of dst.
690+
// That means however, that to the target cells may be assigned values of source
691+
// which are not nearest to the center of dst cells. Usualy probably not a problem
692+
// but we are not precise. The shift is in maximum ... TODO
693+
// b) We could cut the first destination column and left only the second one which is
694+
// completely covered by src. No (significant) stretching is applied in that
695+
// case, but the first column may be rendered as without values event if its center
696+
// is covered by src column. That could result in wrongly rendered (missing) edges
697+
// which could be easily noticed by user
698+
//
699+
// Conclusion: Both are incorrect, but the b) is probably worse because it can be easily
700+
// noticed. Missing data are worse than slightly shifted nearest src cells.
701+
//
702+
703+
if ( xRes > srcXRes )
655704
{
656-
centerRaster = -1. * ( mGeoTransform[3] - center ) / mGeoTransform[5];
705+
if ( mExtent.xMinimum() < myRasterExtent.xMinimum() )
706+
{
707+
srcLeft = qRound(( myRasterExtent.xMinimum() - mExtent.xMinimum() ) / srcXRes );
708+
}
709+
if ( mExtent.xMaximum() > myRasterExtent.xMaximum() )
710+
{
711+
srcRight = qRound(( myRasterExtent.xMaximum() - mExtent.xMinimum() ) / srcXRes ) - 1 ;
712+
}
657713
}
658-
else
714+
if ( yRes > qAbs( srcYRes ) )
659715
{
660-
centerRaster = ( center - mGeoTransform[3] ) / mGeoTransform[5];
716+
if ( mExtent.yMaximum() > myRasterExtent.yMaximum() )
717+
{
718+
srcTop = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMaximum() ) / srcYRes );
719+
}
720+
if ( mExtent.yMinimum() < myRasterExtent.yMinimum() )
721+
{
722+
srcBottom = qRound( -1. * ( mExtent.yMaximum() - myRasterExtent.yMinimum() ) / srcYRes ) - 1 ;
723+
}
661724
}
662-
rasterBottom = static_cast<int>( floor( centerRaster ) );
663-
bottomSpace = myRasterExtent.yMinimum() - ( mGeoTransform[3] + ( rasterBottom + 1 ) * mGeoTransform[5] );
664725

665-
// left
666-
center = myRasterExtent.xMinimum() + xRes / 2;
667-
centerRaster = ( center - mGeoTransform[0] ) / mGeoTransform[1];
668-
rasterLeft = static_cast<int>( floor( centerRaster ) );
669-
leftSpace = myRasterExtent.xMinimum() - ( mGeoTransform[0] + rasterLeft * mGeoTransform[1] );
670-
671-
// right
672-
center = myRasterExtent.xMaximum() - xRes / 2;
673-
centerRaster = ( center - mGeoTransform[0] ) / mGeoTransform[1];
674-
rasterRight = static_cast<int>( floor( centerRaster ) );
675-
rightSpace = ( mGeoTransform[0] + ( rasterRight + 1 ) * mGeoTransform[1] ) - myRasterExtent.xMaximum();
726+
// *** CASE 3: xRes < srcXRes, yRes < srcYRes
727+
// IMHO, we cannot align target grid to raster grid using target grid edges
728+
// and finding the nearest raster grid because it could lead to cell center
729+
// getting outside the right cell when doing resampling, example
730+
// src: | | | |
731+
// dst: | | | |
732+
// Raster width is 30m and it has 3 columns and we want to read xrange 5.1-30
733+
// to 3 columns, the nearest edge for beginning in raster grid is 10.0
734+
// reading cols 1-2, we get raster[1] value in target[0], but the center of
735+
// target[0] is 5.1 + ((30-5.1)/3)/2 = 9.25 so it falls to raster[0]. Is it right?
736+
// => We are looking for such alignment with which the center of first/last cell
737+
// falls to the right raster cell
738+
int xAddPixels = 0;
739+
int leftAddPixels = 0;
740+
if ( xRes < srcXRes )
741+
{
742+
center = myRasterExtent.xMinimum() + xRes / 2;
743+
centerRaster = ( center - mExtent.xMinimum() ) / srcXRes;
744+
srcLeft = static_cast<int>( floor( centerRaster ) );
745+
746+
// Space which has to be added to destination to fit better source
747+
// Warning - TODO: if xRes is much smaller than srcXRes, adding space and pixels to dst may
748+
// result in very large dst grids!!! E.g.
749+
// src: | | |
750+
// dst: | | |
751+
// Solution could be vector rendering of rectangles.
752+
753+
leftSpace = myRasterExtent.xMinimum() - ( mExtent.xMinimum() + srcLeft * srcXRes );
754+
leftSpace = leftSpace > 0 ? leftSpace : 0; // makes only sense for positive
755+
center = myRasterExtent.xMaximum() - xRes / 2;
756+
centerRaster = ( center - mExtent.xMinimum() ) / srcXRes;
757+
srcRight = static_cast<int>( floor( centerRaster ) );
758+
759+
rightSpace = ( mExtent.xMinimum() + ( srcRight + 1 ) * srcXRes ) - myRasterExtent.xMaximum();
760+
rightSpace = rightSpace > 0 ? rightSpace : 0;
761+
QgsDebugMsg( QString( "center = %1 centerRaster = %2 srcRight = %3 rightSpace = %4" ).arg( center ).arg( centerRaster ).arg( srcRight ).arg( rightSpace ) );
762+
763+
double xAdd = leftSpace + rightSpace;
764+
xAddPixels = qRound( xAdd / xRes );
765+
leftAddPixels = qRound( leftSpace / xRes );
766+
}
676767

677-
QgsDebugMsg( QString( "rasterTop = %1 rasterBottom = %2 rasterLeft = %3 rasterRight = %4" ).arg( rasterTop ).arg( rasterBottom ).arg( rasterLeft ).arg( rasterRight ) );
768+
int yAddPixels = 0;
769+
int topAddPixels = 0;
770+
if ( yRes < qAbs( srcYRes ) )
771+
{
772+
center = myRasterExtent.yMaximum() - yRes / 2;
773+
centerRaster = -1. * ( mExtent.yMaximum() - center ) / srcYRes;
774+
srcTop = static_cast<int>( floor( centerRaster ) );
775+
topSpace = ( mExtent.yMaximum() + srcTop * srcYRes ) - myRasterExtent.yMaximum();
776+
topSpace = topSpace > 0 ? topSpace : 0;
777+
center = myRasterExtent.yMinimum() + yRes / 2;
778+
centerRaster = -1. * ( mExtent.yMaximum() - center ) / srcYRes;
779+
780+
srcBottom = static_cast<int>( floor( centerRaster ) );
781+
QgsDebugMsg( QString( "myRasterExtent.yMinimum() = %1 myRasterExtent.yMaximum() = %2" ).arg( myRasterExtent.yMinimum() ).arg( myRasterExtent.yMaximum() ) );
782+
bottomSpace = myRasterExtent.yMinimum() - ( mExtent.yMaximum() + ( srcBottom + 1 ) * srcYRes );
783+
bottomSpace = bottomSpace > 0 ? bottomSpace : 0;
784+
QgsDebugMsg( QString( "center = %1 centerRaster = %2 srcBottom = %3 bottomSpace = %4" ).arg( center ).arg( centerRaster ).arg( srcBottom ).arg( bottomSpace ) );
785+
786+
double yAdd = topSpace + bottomSpace;
787+
yAddPixels = qRound( yAdd / yRes );
788+
topAddPixels = qRound( topSpace / yRes );
789+
}
678790

679-
QgsDebugMsg( QString( "topSpace = %1 bottomSpace = %2 leftSpace = %3 rightSpace = %4" ).arg( topSpace ).arg( bottomSpace ).arg( leftSpace ).arg( rightSpace ) );
791+
srcWidth = srcRight - srcLeft + 1;
792+
srcHeight = srcBottom - srcTop + 1;
680793

681-
int width = right - left + 1;
682-
int height = bottom - top + 1;
794+
QgsDebugMsg( QString( "srcTop = %1 srcBottom = %2 srcLeft = %3 srcRight = %4" ).arg( srcTop ).arg( srcBottom ).arg( srcLeft ).arg( srcRight ) );
683795

684-
int rasterWidth = rasterRight - rasterLeft + 1;
685-
int rasterHeight = rasterBottom - rasterTop + 1;
796+
QgsDebugMsg( QString( "topSpace = %1 bottomSpace = %2 leftSpace = %3 rightSpace = %4" ).arg( topSpace ).arg( bottomSpace ).arg( leftSpace ).arg( rightSpace ) );
686797

687-
QgsDebugMsg( QString( "width = %1 height = %2 rasterWidth = %3 rasterHeight = %4" ).arg( width ).arg( height ).arg( rasterWidth ).arg( rasterHeight ) );
688798

689-
// TODO: what is better floor/ceil, can be negative?
690-
// should be similar
691-
//double xAdd = rasterWidth*rasterXRes - width*xRes;
692-
double xAdd = leftSpace + rightSpace;
693-
int xAddPixels = qRound( xAdd / xRes );
694-
int leftAddPixels = qRound( leftSpace / xRes );
799+
QgsDebugMsg( QString( "width = %1 height = %2 srcWidth = %3 srcHeight = %4" ).arg( width ).arg( height ).arg( srcWidth ).arg( srcHeight ) );
695800

696-
//double leftAdd = rasterWidth*rasterXRes - width*xRes;
697-
double yAdd = topSpace + bottomSpace;
698-
int yAddPixels = qRound( yAdd / yRes );
699-
int topAddPixels = qRound( topSpace / yRes );
700801

701802
QgsDebugMsg( QString( "xAddPixels = %1 yAddPixels = %2 leftAddPixels = %3 topAddPixels = %4" ).arg( xAddPixels ).arg( yAddPixels ).arg( leftAddPixels ).arg( topAddPixels ) );
702-
// Currently only positive allowed, verify if negative has sense and check following use
703-
xAddPixels = xAddPixels > 0 ? xAddPixels : 0;
704-
yAddPixels = yAddPixels > 0 ? yAddPixels : 0;
705-
leftAddPixels = leftAddPixels > 0 ? leftAddPixels : 0;
706-
topAddPixels = topAddPixels > 0 ? topAddPixels : 0;
707803

708804
int totalWidth = width + xAddPixels;
709805
int totalHeight = height + yAddPixels;
710806

711807
QgsDebugMsg( QString( "totalWidth = %1 totalHeight = %2" ).arg( totalWidth ).arg( totalHeight ) );
712808

713-
int size = dataTypeSize( theBandNo ) / 8;
714-
715-
// fill with null values
716-
QByteArray ba = noValueBytes( theBandNo );
717-
char *nodata = ba.data();
718-
char *block = ( char * ) theBlock;
719-
for ( int i = 0; i < thePixelWidth * thePixelHeight; i++ )
720-
{
721-
memcpy( block, nodata, size );
722-
block += size;
723-
}
724809

725810
GDALRasterBandH gdalBand = GDALGetRasterBand( mGdalDataset, theBandNo );
726811
GDALDataType type = ( GDALDataType )mGdalDataType[theBandNo-1];
727812
CPLErrorReset();
728813

729-
// This can be probably used if xAddPixels and yAddPixels are 0 to avoid memcpy
730-
#if 0
731-
// Calc beginnig of data if raster does not start at top
732-
block = ( char * ) theBlock;
733-
if ( top != 0 )
814+
if ( xAddPixels == 0 && yAddPixels == 0 )
734815
{
735-
block += size * thePixelWidth * top;
736-
}
816+
// Calc beginnig of data if raster does not start at top
817+
block = ( char * ) theBlock;
818+
if ( top != 0 )
819+
{
820+
block += dataSize * thePixelWidth * top;
821+
}
737822

738-
// Cal nLineSpace if raster does not cover whole extent
739-
int nLineSpace = size * thePixelWidth;
740-
if ( left != 0 )
741-
{
742-
block += size * left;
823+
// Cal nLineSpace if raster does not cover whole extent
824+
int nLineSpace = dataSize * thePixelWidth;
825+
if ( left != 0 )
826+
{
827+
block += dataSize * left;
828+
}
829+
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
830+
srcLeft, srcTop, srcWidth, srcHeight,
831+
( void * )block,
832+
width, height, type,
833+
0, nLineSpace );
743834
}
744-
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
745-
rasterLeft, rasterTop, rasterWidth, rasterHeight,
746-
( void * )block,
747-
width, height, type,
748-
0, nLineSpace );
749-
#endif
835+
else
836+
{
837+
// Allocate temporary block
838+
void *tmpBlock = malloc( dataSize * totalWidth * totalHeight );
750839

751-
// Allocate temporary block
752-
void *tmpBlock = malloc( size * totalWidth * totalHeight );
840+
CPLErrorReset();
841+
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
842+
srcLeft, srcTop, srcWidth, srcHeight,
843+
( void * )tmpBlock,
844+
totalWidth, totalHeight, type,
845+
0, 0 );
753846

754-
CPLErrorReset();
755-
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
756-
rasterLeft, rasterTop, rasterWidth, rasterHeight,
757-
( void * )tmpBlock,
758-
totalWidth, totalHeight, type,
759-
0, 0 );
847+
if ( err != CPLE_None )
848+
{
849+
QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
850+
QgsDebugMsg( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
851+
free( tmpBlock );
852+
return;
853+
}
760854

761-
if ( err != CPLE_None )
762-
{
763-
QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
764-
QgsDebugMsg( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
765-
free( tmpBlock );
766-
return;
767-
}
855+
for ( int i = 0; i < height; i++ )
856+
{
857+
int r = i + topAddPixels;
858+
char *src = ( char * )tmpBlock + dataSize * r * totalWidth + dataSize * leftAddPixels;
859+
char *dst = ( char * )theBlock + dataSize * ( top + i ) * thePixelWidth + dataSize * ( left );
860+
memcpy( dst, src, dataSize*width );
861+
}
768862

769-
for ( int i = 0; i < height; i++ )
770-
{
771-
int r = i + topAddPixels;
772-
char *src = ( char * )tmpBlock + size * r * totalWidth + size * leftAddPixels;
773-
char *dst = ( char * )theBlock + size * ( top + i ) * thePixelWidth + size * ( left );
774-
memcpy( dst, src, size*width );
863+
free( tmpBlock );
775864
}
776-
777-
free( tmpBlock );
778865
return;
779866
}
780867

0 commit comments

Comments
 (0)
Please sign in to comment.