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
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
