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

[Kay F. Jahnke]

On 18.11.18 13:18, Even Rouault wrote:
> 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.

Yes. Because the prefilter was omitted. The cubic spline interpolation 
is the one method in gdal which needs prefiltering. But b-spline 
interpolation is the best method offered. So I picked it. Then I noticed 
the prefilter is omitted, so I made some noise.

The other interpolators you mention are 'direct' interpolators and 
require no prefiltering.

> (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...)

That's good. I was using

gdalwarp -tr 1 1 -r cubicspline in.tif out.tif, that's where I found the 
difference..

> 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...)

Let's be precise about it. what you link to above is

"Interpolation Revisited" by Philippe Thévenaz, Member, IEEE, Thierry 
Blu, Member, IEEE, and Michael Unser, Fellow, IEEE

from IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 19, NO. 7, JULY 2000

This is the seminal paper which sparked my interest in b-splines some 
years ago. It's really thorough and a lot of it is beyond my 
mathematical horizon, but it sure explains all these things about 
b-splines I've been going on about. P. Thévenaz has also published 
example code for prefiltering and evaluating b-splines here:

http://bigwww.epfl.ch/thevenaz/interpolation/

>> 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 ?

Just an example. I fed my b-spline calculator with a bunch of zeros and 
a single one, and this is what it produced after prefiltering. Another 
example:

-0.464101615138 1.732050807569 -0.464101615138

when you do the prefilter, what you get precisely depends on the *whole* 
signal - that's what is meant by 'indefinite support'. In real life, you 
limit the support to a good handful or two and live with a small error.

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

You may want to stick with CPU implementations for the time being - 
prefiltering for a cubic b-spline is trivial in principle, just dealing 
with the infinite support of the forward-backward recursive filter needs 
some thought, especially when no-data values are involved. Have a look 
at Thévenaz' code I've linked to above, you'll see it's not difficult. 
If you need help with it, ask me.

There is no lack of papers on b-splines. This is very well known territory.

>> 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 ;-)

I *may* consider this.

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

I think here you are - at least partly - mistaken. Let me refer you to

https://www.gdal.org/gdal_vrttut.html

" ... For now kernel is not applied to sub-sampled or over-sampled data."

This makes sense: you can only apply a convolution to discrete samples. 
It's not an interpolator. So in order to subsample, you have to first 
use the VRT file to establish the convolution kernel, then save to a new 
file to have it applied, and then subsample the result with an 
interpolator. This makes the kernel-filtered source unsuitable for 
subsampling - unless it's subsampling by integer-valued factors 
('decimation'), where it could be applied on-the-fly, a special case 
which gdal does not seem to exploit.

Before I call it a day - let's not forget my core issue, which is the 
necessity for low-pass-filtering before subsampling. The -r cubicspline 
is a side issue which I noticed on top of the missing low-pass, as are 
my observations about -r average and 'area interpolation'. You can avoid 
these two issues by picking other interpolators in gdal, but you can't 
avoid the low-pass filter before subsampling.

Kay