On 09/19/2014 01:50 PM, Peter Hopfgartner wrote:
On 09/18/2014 10:30 PM, Even Rouault wrote:
Le jeudi 18 septembre 2014 22:23:06, Peter Hopfgartner a écrit :
Salut Evan,
thanks for you recipe. Unfortunately, what I am really interested in
is to
use this from MapServer. Could you give me a pointer in the source
where
the resampling takes place?
In gcore/rasterio.cpp : GDALRasterBand::IRasterIO(). I must warn you
: quite
tricky code
But perhaps you don't need code change. In MapServer did you try to
define
PROCESSING "RESAMPLE=BILINEAR" ?
Thanks,
Peter
Am 18.09.2014 um 22:11 schrieb Even Rouault
<even.roua...@spatialys.com>:
Le jeudi 18 septembre 2014 17:40:11, Peter Hopfgartner a écrit :
Hi, everybody
we are using the GDAL WMS driver for getting maps from a public WMTS
service, http://www.basemap.at. Anyway, the resulting image does not
look smooth. In particular, I would guess that the the image was
scaled
with something like a nearest neighbor algorithm.
So my question is simply, if I can switch the interpolation method in
some place.
The XML file, taken from
http://gisforge.wordpress.com/2012/11/06/basemap-at-wmts-layer-in-qgis/,
is:
<GDAL_WMS>
<Service name="TMS">
<ServerUrl>http://maps.wien.gv.at/basemap/geolandbasemap/normal/google38
57/ ${z}/${y}/${x}.jpeg</ServerUrl> </Service>
<DataWindow>
<UpperLeftX>-20037508.34</UpperLeftX>
<UpperLeftY>20037508.34</UpperLeftY>
<LowerRightX>20037508.34</LowerRightX>
<LowerRightY>-20037508.34</LowerRightY>
<TileLevel>18</TileLevel>
<TileCountX>1</TileCountX>
<TileCountY>1</TileCountY>
<YOrigin>top</YOrigin>
</DataWindow>
<Projection>EPSG:3857</Projection>
<BlockSizeX>256</BlockSizeX>
<BlockSizeY>256</BlockSizeY>
<BandsCount>3</BandsCount>
<Cache />
</GDAL_WMS>
The command is:
gdal_translate -srcwin 35500000 23500000 100000 100000 -of JPEG
-outsize
250 250 basemap_at.xml basemap_at.jpg
Peter,
Resampling done by gdal_translate is nearest neighbour, so not always
that great.
What you can do is select an outsize that is compatible with the
source
overviews. In TMS, the ratio between sucessive overviews work with
factor
of 2, so you should select an outsize of 10000 / 2^N. The closest is
10000 / 256 = 390.625 -> 391
So
gdal_translate -srcwin 35500000 23500000 100000 100000 -outsize
391 391 basemap_at.xml tmp.tif
and then:
gdalwarp tmp.tif tmp2.tif -r lanczos -ts 250 250 (you can try
different
resampling kernel)
gdal_translate tmp2.tif basemap_at.jpg -of JPEG
Best regards,
Even
I've tried to simulate the result of different algorithms to further
investigate the idea of having a configurable resampling step. The
original tile is 256 x 256 pixels and I've rescaled it to rescale it
to 200x200 using Gimp.
The original tile was taken from
http://maps.wien.gv.at/basemap/geolandbasemap/normal/google3857/8/89/137.jpg,
a copy is stored at
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.jpg.
Results
no interpolation, (nearest neighbor?):
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-noint.png
bilinear:
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-linear.png
bicubic:
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-cubic.png
lanczos:
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-lanczos.png
I'd say that the quality of the resampling respects the obove order,
with lanczos being the best.
Regards,
Peter
I finally took the time to have a better look at this.
The attached patch is only for our special case, where data arrives in 8
bit bands and is certainly far from being decent (I've written my last
C++ code approx. 15 years ago). Anyway, the image quality seems indeed
to improve sensibly.
I used the current SVN version as a base for the development.
Applying the patch I get from
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/basemap.2000x2000-quasiorig.jpg
to this
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/basemap.2000x2000.bilinear.jpg
.
Anyway, the execution time approximatly doubles (on my 2011 notebook it
passes from ~ 1.8 s to ~ 3.4 s), which in our case is not really an
issue, but probaly for almost everybody else.
The command line is, starting from the compilation dir:
apps/gdal_translate -srcwin 35555000 23538000 2500 2500 -of JPEG
-outsize 2000 2000 basemap.xml basemap.2000x2000.bilinear.jpg .
Any comment is highly appreciated.
Regards and a nice weekend,
Peter
--
Peter Hopfgartner
R3 GIS
web : http://www.r3-gis.com
Index: gcore/rasterio.cpp
===================================================================
--- gcore/rasterio.cpp (revision 27739)
+++ gcore/rasterio.cpp (working copy)
@@ -70,7 +70,7 @@
GByte *pabySrcBlock = NULL;
GDALRasterBlock *poBlock = NULL;
int nLBlockX=-1, nLBlockY=-1, iBufYOff, iBufXOff, iSrcY;
-
+
if( eRWFlag == GF_Write && eFlushBlockErr != CE_None )
{
CPLError(eFlushBlockErr, CPLE_AppDefined,
@@ -368,7 +368,7 @@
for( iDstX = nXOff; iDstX < nXOff + nXSize; iDstX ++)
{
- iBufXOff = (int)((iDstX - nXOff) / dfSrcXInc);
+ iBufXOff = (int)((iDstX - nXOff) / dfSrcXInc);
iBufOffset = (size_t)iBufYOff * nLineSpace + iBufXOff * nPixelSpace;
/* -------------------------------------------------------------------- */
@@ -448,10 +448,10 @@
dfSrcY = (iBufYOff+0.5) * dfSrcYInc + nYOff;
dfSrcX = 0.5 * dfSrcXInc + nXOff;
+ double dFracY = fmod(dfSrcY, 1.0);
iSrcY = (int) dfSrcY;
iBufOffset = (size_t)iBufYOff * nLineSpace;
-
if( iSrcY >= nLimitBlockY )
{
nLBlockY = iSrcY / nBlockYSize;
@@ -462,12 +462,22 @@
nStartBlockX = -nBlockXSize; /* make sure a new block is loaded */
size_t iSrcOffsetCst = (iSrcY - nLBlockY*nBlockYSize) * (size_t)nBlockXSize;
+ size_t iSrcOffsetCstTop;
+
+ size_t iSrcYTop;
+ if ( iSrcY < nLimitBlockY - 1 ) {
+ iSrcYTop = iSrcY + 1;
+ iSrcOffsetCstTop = (iSrcYTop - nLBlockY*nBlockYSize) * (size_t)nBlockXSize;
+ } else {
+ iSrcYTop = iSrcY;
+ iSrcOffsetCstTop = iSrcOffsetCst;
+ }
for( iBufXOff = 0; iBufXOff < nBufXSize; iBufXOff++, dfSrcX += dfSrcXInc )
{
iSrcX = (int) dfSrcX;
+ double dFracX = fmod(dfSrcX, 1.0);
int nDiffX = iSrcX - nStartBlockX;
-
/* -------------------------------------------------------------------- */
/* Ensure we have the appropriate block loaded. */
/* -------------------------------------------------------------------- */
@@ -494,15 +504,37 @@
return CE_Failure;
}
}
+ int nDiffXRight;
+ if (nDiffXRight < nBlockXSize - 1) {
+ nDiffXRight = nDiffX + 1;
+ } else {
+ nDiffXRight = nDiffX;
+ }
/* -------------------------------------------------------------------- */
/* Copy over this pixel of data. */
/* -------------------------------------------------------------------- */
iSrcOffset = ((size_t)nDiffX + iSrcOffsetCst)*nBandDataSize;
-
if( bByteCopy )
{
- ((GByte *) pData)[iBufOffset] = pabySrcBlock[iSrcOffset];
+ size_t iSrcOffsetBR = ((size_t)nDiffXRight + iSrcOffsetCst)*nBandDataSize;
+ size_t iSrcOffsetTL = ((size_t)nDiffX + iSrcOffsetCstTop)*nBandDataSize;
+ size_t iSrcOffsetTR = ((size_t)nDiffXRight + iSrcOffsetCstTop)*nBandDataSize;
+
+ unsigned char ucVal;
+ if (eDataType == GDT_Byte) {
+ unsigned char v00, v01, v10, v11;
+ v00 = pabySrcBlock[iSrcOffset];
+ v01 = pabySrcBlock[iSrcOffsetBR];
+ v10 = pabySrcBlock[iSrcOffsetTL];
+ v11 = pabySrcBlock[iSrcOffsetTR];
+
+ double valY0 = double(v00) * (1 - dFracX) + double(v01) * dFracX;
+ double valY1 = double(v10) * (1 - dFracX) + double(v11) * dFracX;
+ double dVal = valY0 * (1 - dFracY) + valY1 * dFracY;
+ ucVal = (unsigned char)round(dVal);
+ }
+ ((GByte *) pData)[iBufOffset] = ucVal;
}
else if( eDataType == eBufType )
{
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
