I've read some articles which introduce fast third-order interpolation using GL_LINEAR.
Because [1] contains lots of errata, I recommend reading [2] if you want to catch the formalism.
Both of them mention the restriction of this method. For filtered texture with GL_LINEAR, next relation holds only if 0 <= b/(a+b) <= 1
a*f(i, j) + b*f(i+1, j) = F(i+b/(a+b), j)
where f is original image data and F is linearly interpolated texture by OpenGL.
Here's the problem. [1] mentions that this method also can be applied to Catmull-Rom bicubic.
the method can also be adapted to interpolating filters such as Catmull-Rom splines
However, it's obvious that with Catmull-Rom weighting function which contains negative parts, the condition(0 <= b/(a+b) <= 1) cannot be fullfiled. In fact, I have tried to implement Catmull-Rom with same logic, it produces just blurry images.
Is there a special way to apply the method in [1] and [2] to Catmull-Rom interpolation? Or Do I have to fetch all 16 texels for Catmull-Rom?
I think that the paper by Sigg and Hadwiger is correct.
Catmull-Rom polynomial p(t) can be written as
p(t) = 0.5 [ w0(t) * p0 + w1(t) * p1 + w2(t) * p2 + w3(t) * p3]
where
w0(t) = -t + 2*t^2 - t^3;
w1(t) = 2 - 5*t^2 + 3*t^3;
w2(t) = t + 4*t^2 - 3*t^3;
w3(t) = -t^2 + t^3;
and p0, p1, p2, p3 are the sampled function values.
For the B-spline, you group together the first two terms, namely w0(t) * p0 and w1(t) * p1 and the second two terms, namely w2(t) * p2 and w3(t) * p3. For the Catmull-Rom spline, you group together w0(t) * p0 and w3(t) * p3 and w1(t) * p1 and w2(t) * p2. It is easy to verify, even by a Matlab plot, that the condition b/(a+b) is satisfied for this choice. So, the idea can be used also in Catmull-Rom interpolation. However, for the Catmull-Rom case, it seems to me that there is no possible extension to the last step, since neither (w0(t) + w3(t))/(w0(t) + w1(t) + w2(t) + w3(t)) nor (w1(t) + w2(t))/(w0(t) + w1(t) + w2(t) + w3(t)) meet the condition you specified.