Bicubic B-Spline interpolation (the Bourke algorithm)

[mschuster]

I have been playing with StarAlignment's various interpolation options and have concluded that the Bourke algorithm used by PI (according to the documentation) is not a correct implementation of Bicubic B-Spline interpolation.

The Bourke algorithm uses the Bicubic B-Spline basis function but does not weight instances of this basis properly. The Bicubic B-spline basis has no negative lobes, so using this basis directly with pixel weights as Bourke does basically amounts to a convolution with a poor approximation to a Gaussian. This accounts for the blurry results.

The correct basis weights are not pixel values themselves but rather a transform of them. The transform can be implemented in linear time either by a banded matrix algorithm or by a digital filtering algorithm. Basically in pixel space the support of the B-spline interpolator is infinite (surprisingly) and the transform provides an efficient implementation.

I implemented the digitial filtering algorithm using the references below. It gives sharp results, much like those of the higher order Lanczos interpolators. Unfortunately on my strongly undersampled images there too much ringing around stars with FWHM approaching 1 pixel. In the smoother nebula areas however the results are very good.

Here are some references:

Unser et al., "Fast B-Spline Transforms for Continuous Image Representation and Interpolation", IEEE Transactions on Pattern Analysis and Machine Intelligence, 13(3):277-285, 1991 March.

Condat et al., "Beyond Interpolation: Optimal Reconstruction by Quasi-Interpolation", International Conference on Image Processing 2005, pp. 33-36.

Thanks,
Mike

[cs_pixinsight]

Mike, Wow!  Your determination in tracking this down is amazing.  This finding will definitely help PI out for Image Integration issues with Lanczos dark ringing.  Hopefully Juan will concur with your findings and be able to implement an algorithm change.

Craig

[mschuster]

Craig,

Just to be clearer and trying to avoid confusion, I think most of us are using either "Lanczos" or "Bicubic spline" or "Auto", but not "Bicubic B-spline".

The names are similar so confusion is likely. It is "Bicubic B-spline", the one with the "B-", that I think there is a problem.

The first three are working fine I believe (although there was a bug report a while ago, but its probably fixed now).

The dark ringing is a different issue probably. And actually the "Bicubic B-spline" fix I am suggesting, which should increase sharpness, will make dark ringing more likely.

Mike

[ruediger]

Mike, thanks for this post!

In a recent project I used Bicubic B-Spline the first time and noticed the blurry stacked result, but had absolutely no idea what's causing this, because I also changed other parameters in the preprocessing pipeline as well. So now it's time to register 301 frames again... this time with Bicubic spline interpolation. I visually compared some Bicubic spline registered frames with Bicubic B-Spline frames and the difference in sharpness is directly seen.

Rüdiger

[Juan Conejero]

Hi Mike,

I wouldn't say that PI's implementation of Bourke's interpolation is a 'bug'. It is just an implementation of a well-known ('classical', so to say) interpolation algorithm 'as is'. Obviously, it is now rarely used due to its very poor performance, and it's part of PI/PCL mostly for completeness and backward-compatibility reasons (it was also available in PI LE).

In the image interpolation problem we always have to find a delicate balance between detail preservation and prevention of artifacts. This is particularly problematic for interpolation of linear data, due to wild variations between neighbor pixels. A positive interpolation function (without negative lobes) behaves like a low-pass filter, and hence does not generate ringing. However, it does generate aliasing and tends to blur the data. A function with negative lobes, such as cubic splines or the Lanczos functions, allows us to control the exact dose of high-pass filtering to improve detail preservation and reduce aliasing, but at the cost of some ringing. In our implementations we use clamping methods to allow for a finer control of the ringing versus aliasing problem. The perfect interpolation algorithm does not exist; the user has to learn to analyze the data and find the best compromise for each case. We only can try to provide a comprehensive set of the best possible implementations.

I encourage you to explore and implement new interpolation algorithms. As you know, PCL is an open-source library since version 2.0, so this is an excellent opportunity to contribute. I already knew the references you have included, and IMO the quasi-interpolators are particularly interesting.

[mschuster]

Juan FYI: I played with the quasi-interpolators. They look good but all ring on my subs (quasi-linear, quasi-quadratic and quasi-cubic). My attempts at clamping don't work well I think because the ringing extends several pixels outward from the pair of adjacent pixels that exceed the clamping threshold.
Mike