Third-Order Texture Filtering using Linear Interpolation

A straightfoward implementation of bi- or tricubic texture filtering requires lots of texture lookups. As texture lookups are slow on current graphics hardware, reducing their number greatly improves the performance of high-order filtering kernels. In contrast, linear interpolation is relatively cheap as it is supported by the hardware. This fact has been exploited by Sigg and Hadwiger (2005) to make cubic interpolation much faster by performing a series of linear texture lookups. In this article, an in-depth explanation of their technique is given, starting with some general information on linear interpolation and approximation using cubic B-splines.

At first, we will focus on the one dimensional case. Suppose we have two data points $f_i$ and $f_{i+1}$ and we want to interpolate the values between these points. Using linear interpolation, this results in the following equation:

$f(t) = w_0(t) f_i + w_1(t) f_{i+1} = (1-t) f_i + t f_{i+1}$ ,

where $0 \le t \le 1$ . The functions $w_0(t)=1-t$ and $w_1(t)=t$ determine the filtering weight for each point. The values are interpolated such that $f(t)$ is $f_i$ at $t=0$ and $f_{i+1}$ at $t=1$ . The following images show the liner basis functions (left) and the result of interpolating some points (right):

The formula is automatically evaluated on the GPU if linear filtering is activated in OpenGL and, for example, texture is used in the fragment shader to fetch a texel from texture memory.

If we want to use cubic filtering, we have to use the following formula:

$f(t) = w_0(t) f_{i-1} + w_1(t) f_i + w_2(t) f_{i+1} + w_3(t) f_{i+2}$ ,

where $w_j(t)$ are now third-order polynomials. We will choose weightings based on uniform cubic B-splines, so the weighting functions become

$w_0(t) = \frac{1}{6} (-t^3 + 3t^2 - t + 1) \\ w_1(t) = \frac{1}{6} (3t^3 - 6t^2 + 4) \\ w_2(t) = \frac{1}{6} (-3t^3 + 3t^2 + 3t + 1) \\ w_3(t) = \frac{1}{6} t^3$ .

Although four points $f_{i-1} ... f_{i+2}$ are used, the function constructs only points between $f_i$ and $f_{i+1}$ for $0 \le t \le 1$ . This range is called a segment of the curve. However, rather than interpolating the points $f_{i}$ and $f_{i+1}$ as above, $f(t)$ only approximates the values now. That means that at $t=0$ and $t=1$ , $f(t)$ does not pass through the data points $f_i$ and $f_{i+1}$ . Nevertheless, the resulting function is geometrically continuous, as the approximated value at $t=1$ of a segment is the same as the approximated value at $t=0$ of the next segment. This is shown on the following image on the right (left: cubic B-spline basis functions). Note that the curve does not pass through the blue dots and that the first and last segment is missing.

As we can see from the equation for cubic interpolation, four texture lookups are necessary. This becomes even worse in two and three dimensions, where 16 and 64 lookups are required, respectively. As mentioned in the introduction, it is possible to replace a cubic interpolation by a sequence of linear interpolations. As a result, in 1D the 4 lookups can be reduced to 2 linear lookups. The 16 lookups in 2D can be replaced by 4 linear lookups, the 64 lookups in 3D by 8. In general, $4^n$ lookups are replaced by $2^n$ linear lookups.

Now, let’s focus on how the optimization works. Actually, it is pretty simple: just write the cubic interpolation as the sum of two linear interpolations:

$f(t) = \underbrace{w_0(t) f_{i-1} + w_1(t) f_i}_{\text{First linear interpolation}} + \underbrace{w_2(t) f_{i+1} + w_3(t) f_{i+2}}_{\text{Second linear interpolation}}$ .

For this purpose, let’s first find a notation for the operation that the graphics card executes when performing a linear texture lookup:

$f(t) = (1-t) f_i + t f_{i+1} = f_{i+t}$ .

That shows that when we want to interpolate between two texels, we just extract a constant from texture memory at $f_{i+t}$ . However, that is not sufficient if we want to write weighted sums $a f_{i} + b f_{i+1}$ with two weights $a$ and $b$ as linear interpolation. So, we have to scale $a$ and $b$ , resulting in the new variables $a'$ and $b'$ with $a'+b'=1$ :

$a' = \frac{a}{a+b} \\ b' = \frac{b}{a+b}$ .

Then, we write the sum as follows:

$(a+b) (a' f_{i} + b' f_{i+1}) = (a+b) ((1-b') f_{i} + b' f_{i+1})$ .

Notice that this looks exactly like an ordinary linear interpolation, scaled by the factor $a+b$ . Using the above notation, the weighted sum becomes

$a f_{i} + b f_{i+1} = (a+b) f_{i+\frac{b}{a+b}}$ .

There is a small restriction: it only works if $0 \le \frac{b}{a+b} \le 1$ . Fortunately, that is satisfied when using cubic B-splines.

Finally, third-order texture filtering is obtained by the following two linear interpolations:

$f(t) = (w_0(t) + w_1(t)) f_{i-1+\frac{w_1(t)}{w_0(t) + w_1(t)}} \\ ~~~~~~ + (w_2(t) + w_3(t)) f_{i+1+\frac{w_3(t)}{w_2(t) + w_3(t)}}$ .

In conclusion, this post has presented an efficient technique for third-order texture filtering in 1D based on hardware-accelerated linear interpolation. In the next post, this will be extended to two and three dimensions, where the performance gain is much more significant.

References:

Sigg and Hadwiger (2005): Fast Third-Order Texture Filtering. In GPU Gems 2: Programming Techniques for High-Performance Graphics and General-Purpose Computation, edited by Matt Pharr, pp. 313–329.

0xef