Lanczos Interpolation

The Lanczos filter

The Oxford English Dictionary defines interpolation as "the act of adding a value into a series by calculating it from surrounding known values".

For example, we have a set of values:

\[ \begin{align} &S = (2, 0, 1.5, 1)\\\\ &s_0 = 2,\quad s_1 = 0,\quad s_2 = 1.5,\quad s_3 = 1 \end{align} \]

and we add a value at \(x=1.4\), between \(s_1\) and \(s_2\):

\[ \begin{align} &S = (2, 0, ?, 1.5, 1)\\\\ &s_1 = 0,\quad s_{1.4} =\:?,\quad s_2 = 1.5 \end{align} \]

The interpolation value depends on the interpolation method. For example, with linear interpolation the value is \(0.6\):

To interpolate a value with Lanczos interpolation, we do a weighted sum of the surrounding values:

\[ S(x) = \displaystyle \sum_{i=\lfloor x \rfloor - a + 1}^{\lfloor x \rfloor + a} s_i L(x - i) \]

Where \(\lfloor x \rfloor\) is the floor function; \(a\) is a possitive integer (2 throughout this example); \(L(x)\) is the weight function (the Lanczos filter):

\[ L(x) = \begin{cases} 1 & \quad \text{if } x = 0\\ \frac{a\:sin(\pi x) sin(\pi x / a)}{\pi^2 x^2} & \quad \text{if } -a \leq x < a\\ 0 & \quad \text{otherwise.} \end{cases} \]

This function in Python would be like:

from numpy import pi, sin

a = 2

def L(x):
    if x == 0: return 1
    elif -a <= x < a:
        return a*sin(pi*x)*sin(pi*x/a) / (pi**2 * x**2)
    else: return 0

For example, the value at \(1.4\) gives us:

\[ \begin{align} &S(1.4) = \displaystyle \sum_{i=0}^{3} s_i L(1.4 - i)\\ &S(1.4) = s_0 L(1.4) + s_1 L(0.4) + s_2 L(-0.6) + s_3 L(-1.6)\\ &S(1.4) \approx 0.4463 \end{align} \]

Python code:

from numpy import pi, sin, floor

a = 2

def L(x):
    if x == 0: return 1
    elif -a <= x < a:
        return a*sin(pi*x)*sin(pi*x/a) / (pi**2 * x**2)
    else: return 0


# ys: y-coordinates
# x:  any value between 0 and len(ys) - 1
def lanczos_query(ys, x):
    n = ys.shape[0]

    if x < 0 or x > n-1:
        raise ValueError('value x must be between 0 and n-1')

    y = 0.0

    for i in range(int(floor(x)-a+1), int(floor(x)+a+1)):
        j = i

        if i < 0: j = 0
        elif i > n-1: j = n-1

        y += ys[j] * L(x - i)

    return y

Watch the video: