[gdal-dev] downsampling geotiff with a low-pass filter

[Even Rouault]

Hi,

> Let' stick with the cubic b-spline for now, and do a sample evaluation
> at an integral position (so, at a sample loaction). Let the clean signal
> be ... 0 , 1 , 0 ... and we'll evaluate right where the '1' is.
> 
> The b-spline reconstruction kernel (or, the set of weights at integral
> positions) is 1/6 2/3 1/6. If we convolve, we get
> 
> 0 * 1/6 + 1 * 2/3 + 0*1/6 = 2/3

OK, indeed experimenting with "gdalwarp -r [method] in.tif out.tif" on in.tif 
having squared pixels, so the evaluation is done exactly at sample locations, 
I found that cubic and lanczos yield to identical result as the input, and 
indeed if you look at the weights, they are 1 at x = 0, and zero at other 
integral position. Wheras cubicspline indeed does not have thir interpolation 
criterion.
(note if you do it with gdal_translate -r without resizing, GDAL internals 
will see that no resampling is needed, and completely by-pass all the 
resampling code, so you'd get apparent perfect reconstruction...)

Found this paper http://bigwww.epfl.ch/publications/thevenaz0002.pdf 
confirming what you say :
"""No B-spline of degree n > 1 benefits from the property of being 
interpolating; [...] Unfortunately, some authors have failed to observe this 
rule, which led them to claim that B-splines typically blur data.  There-
fore,  such  claims  are  misleading,  even  though  they  can
be repeatedly found in the published literature"""
I'm not sure how widely cubicspline is used by GDAL users, but I suspect that 
some people avoid it because of this blurring side-effect (or conversely, 
maybe use it because they find it useful for their use case. who knows...)

> 
> You can see that the interpolation criterion is *not* fulfilled, because
> we'd expect to see '1' as the result. Why is this? because we haven't
> prefiltered. I'll not go into prefiltering now, but a *prefiltered*
> signal at this position is
> 
>   -0.464285714 1.732142857 -0.464285714

Where does those magic values come from ? 

Found this paper speaking about prefiltering for cubic b-spline interpolation:
http://dannyruijters.nl/docs/cudaPrefilter3.pdf

> I'd like to help, but rather by explaining things. I'd not really like 
> to get involved in coding, as I have a lot of work on my hands already.

Documentation contributions are also accepted ;-)

Relevant file:
https://github.com/OSGeo/gdal/blob/master/gdal/apps/gdal_utilities.dox
https://github.com/OSGeo/gdal/blob/master/gdal/apps/gdalwarp_bin.cpp
https://github.com/OSGeo/gdal/blob/master/gdal/gcore/gdal.h#L128
https://github.com/OSGeo/gdal/blob/master/gdal/alg/gdalwarper.h#L48

Although we might want to centralize a bit detailed information about the 
properties of the algorithm if we go to detail them.


> Jukka's hint to use kernel-filtered sources in a VRT file does do the
> trick and shows the way to go, but requires an intermediate file of the
> same size as the initial input. 

No. If you use the VRT directly, it is just a small XML file describing the 
processing you want to apply on source(s). You can use it directly as the 
input of gdal_translate -r, and the VRT specific computations will be done on-
the-fly.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com