[gdal-dev] Interpolation when using the WMS format

[Peter Hopfgartner]

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.rouault at 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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: rasterio.bilinear.patch
Type: text/x-patch
Size: 4284 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140926/8a648293/attachment.bin>