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