Description of bicubic interpolation as implemented in octave's interp2.m.

bicubic.tex

\documentclass{scrarticle}

\usepackage{listings}
\lstset{%
  language=Matlab,
  columns=flexible,
%   flexiblecolumns=true,
  basicstyle=\ttfamily\small,
  commentstyle=\color{green!80!black}\itshape,
  keepspaces=true,
  breaklines=false,
}

% \usepackage{polyglossia} % replaces babel
% \setdefaultlanguage{english}

\usepackage{tikz}
\usepackage{amsmath}
\usepackage{mathtools}
\usepackage{hyperref}

\title{Bicubic interpolation algorithm}
\author{Christof Kaufmann}

% one infinite page, https://tex.stackexchange.com/a/19241/115883
\usepackage[paperheight=\maxdimen, textwidth=17.5cm]{geometry} % width corresponds to DIV12
\usepackage[active,tightpage]{preview}
\renewcommand{\PreviewBorder}{1in}

% preserve parindent, https://tex.stackexchange.com/a/98214
\usepackage{etoolbox}
\edef\keptparindent{\the\parindent}
\patchcmd{\preview}
  {\ignorespaces} %%% \preview ends with \ignorespaces
  {\parindent\keptparindent\ignorespaces}
  {}{}
\begin{document}
\begin{preview}

\maketitle



\section{Inner points}

Bicubic interpolation can be separated into $x$ axis and $y$ axis. For each axis a convolutional algorithm can be applied. The $xy$ grid must be rectilinear, i.\,e.\ $x_{n+1} - x_n = \delta x$ and $y_{n+1} - y_n = \delta y$ must be the same for all $n$, but $\delta x$ and $\delta y$ may differ, and may even be negative. The cubic kernel function with $a = -\frac 1 2$ is according to \href{https://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm}{Wikipedia}:

\noindent
\begin{minipage}{0.38\linewidth}
  \begin{tikzpicture}[x=9mm, y=9mm]
    \draw [->, thick] (-3, 0) -- (3, 0) node [below] {$h$};
    \draw [->, thick] (0, -0.15) -- (0, 1.25) node [right, blue] {$W(h)$};
    \draw [gray] (0.1, 1) -- (-0.1, 1) node [left] {$1$};
    \foreach \x in {-2, ..., 2}
      \draw [gray] (\x, 0.1) -- (\x, -0.1) node [below] {$\x$};
    \draw [blue] (-3, 0) -- (-2, 0)
                         plot [domain=-2:-1] ({\x}, {0.5 * \x*\x*\x + 2.5 * \x*\x + 4 * \x + 2})
                         plot [domain=-1:1] ({\x}, {1.5 * abs(\x*\x*\x) - 2.5 * \x*\x + 1})
                         plot [domain=1:2] ({\x}, {-0.5 * \x*\x*\x + 2.5 * \x*\x - 4 * \x + 2})
                         -- (3, 0);
  \end{tikzpicture}
\end{minipage}%
\begin{minipage}{0.62\linewidth}
  \begin{equation}
    W(h) = 
    \begin{cases}
      \frac 3 2 \, |h|^3 - \frac 5 2 \, |h|^2 + 1 & \text{if } |h| \le 1\\
      -\frac 1 2 \, |h|^3 + \frac 5 2 \, |h|^2 - 4 \, |h| + 2 & \text{if } 1 < |h| \le 2\\
      0 &\text{otherwise}
    \end{cases}
  \end{equation}
\end{minipage}

\noindent
The corresponding code is straightforward:
\begin{lstlisting}[gobble=2]
  function w = cubic(h)
    absh = abs(h);
    absh01 = absh <= 1;            % for |h| <= 1
    absh12 = absh <= 2 & ~absh01;  % for 1 < |h| <= 2
    w = (1.5 * absh.^3 - 2.5 * absh.^2 + 1)             .* absh01 + ...
        (-0.5 * absh.^3 + 2.5 * absh.^2 - 4 * absh + 2) .* absh12;
  end
\end{lstlisting}

\noindent
With this formula it is easy to interpolate a single location using two neighbors around:
\begin{equation}
  y_i = W(-1-h) \, y_0 + W(0-h) \, y_1 + W(1-h) \, y_2 + W(2-h) \, y_3
\end{equation}
Let's illustrate this with a simple 1D example:\\
\begin{minipage}[t]{0.62\linewidth}
  \begin{lstlisting}[gobble=4, belowskip=0pt, aboveskip=0pt]
    x = 1:4;
    y = [2, 1, 0.5, 1.5];
    xi = 2.25;
    
    % for indexes startig at 1:
    idx = floor(xi); 
    h = xi - idx;    
    % or more general:
    idx = lookup (x, xi);
    h = (xi - x(idx)) / abs(x(1) - x(2));

    % as single statement:
    yi = cubic(-1-h)*y(-1+idx) + cubic(0-h)*y(0+idx) ...
       + cubic(1-h)*y(1+idx) + cubic(2-h)*y(2+idx);
    % or as for loop:
    yi = 0;
    for shift = -1:2
      yi += cubic(shift-h) * y(shift+idx);
    endfor
    % or vectorized:
    shifts = -1:2;
    yi = sum(cubic(shifts - h) .* y(shifts + idx));
  \end{lstlisting}
%   yi = interp1(x, y, xi, 'cubic');
\end{minipage}%
\hspace{-0.23\linewidth}%
\begin{minipage}[t]{0.6\linewidth}
  \vspace{0em}
  \begin{tikzpicture}[scale=2]
    \draw [<->, thick] (0, 1.75) |- (4.25, 0);
    \foreach \x/\y [count=\i from 0] in {1/0.25, 2/1, 3/0.5, 4/1.5} {
      \fill (\x, \y) circle [radius=1pt];
      \draw (\x, 0.1) -- (\x, -0.1) node [below] {$x_\i$};
      \draw (0.1, \y) -- (-0.1, \y) node [left] {$y_\i$};
    }
    \draw [blue] (2.25, 103/128) circle [radius=1.5pt];
    \draw [blue] (2.25, 0.1) -- ++(down:0.2) node [below] {$x_i$};
    \draw [blue] (0.1, 103/128) -- ++(left:0.2) node [left] {$y_i$};
    
    \draw [gray, |<->|] (2, 0.1) -- (2.25, 0.1) node [above, midway] {$h$};
  \end{tikzpicture}
\end{minipage}
\bigskip

\noindent
This can be extended to 2D easily:\\
\begin{minipage}{0.28\linewidth}
  \centering
  \begin{tikzpicture}[yscale=-1, x=12mm, y=12mm, font=\ttfamily]
    \foreach \y in {1, ..., 5} {
      \foreach \x in {1, ..., 5} {
        \fill (\x, \y) circle [radius=1.5pt];
      }
      \foreach \x in {2.25, 3.75} {
        \draw [red] (\x, \y) circle [radius=2pt];
      }
    }
    \foreach \y in {2.75, 3.75} {
      \foreach \x in {2.25, 3.75} {
        \draw [blue] (\x, \y) circle [radius=2pt];
      }
    }
    \node at (1.5, 1.5) {A};
    \node [red] at (3, 1.5) {B};
    \node [blue] at (3, 3.5) {C};
  \end{tikzpicture}
\end{minipage}%
\hfill%
\begin{minipage}{0.67\linewidth}
  \begin{lstlisting}[gobble=4]
    A = magic(5)
    xi = [2.25, 3.75];
    yi = [2.75; 3.75];

    idx = floor(xi);
    h = xi - idx;
    B = cubic(-1-h) .* A(:, -1+idx) + cubic(0-h) .* A(:, 0+idx) ...
      + cubic(1-h) .* A(:, 1+idx) + cubic(2-h) .* A(:, 2+idx)

    idx = floor(yi);
    h = yi - idx;
    C = cubic(-1-h) .* B(-1+idx, :) + cubic(0-h) .* B(0+idx, :) ...
      + cubic(1-h) .* B(1+idx, :) + cubic(2-h) .* B(2+idx, :);
  \end{lstlisting}
%   C = interp2(A, xi, yi, "cubic");
\end{minipage}



\section{Border points}

At the border, there are too less neighbors. To apply the algorithm, the matrix can be padded. The simplest padding method would be to add 0s, which is not done in Octave. Still, different approaches are used in different functions.

In \lstinline|imresize| even \lstinline|xi| less than $1$ (but always greater than $0.5$) and thus \lstinline|idx| equal $0$ can occur. Similarly, too large indexes can occur. These points lack two neighbors, but interpolation is still performed. For that, \lstinline|imresize| uses symmetrical padding.

In \lstinline|interp2| with bicubic interpolation any \lstinline|xi| less than $1$ or larger than the number of columns of the matrix \lstinline|z| will yield a \lstinline|NaN| or constant extrapolation value \lstinline|extrap|. So there can only occur situations where one neighbor is missing. In \lstinline|interp2| the missing neighbor is reconstructed using quadratical extrapolation. The cubic interpolation will then degenerate to a quadratic interpolation. The quadratic Lagrange interpolation formula is:
\begin{equation}
  q(x) = \frac{(x - x_2) \, (x - x_3)}{(x_1 - x_2) \, (x_1 - x_3)} \, y_1
       + \frac{(x - x_1) \, (x - x_3)}{(x_2 - x_1) \, (x_2 - x_3)} \, y_2
       + \frac{(x - x_1) \, (x - x_2)}{(x_3 - x_1) \, (x_3 - x_2)} \, y_3
  \label{eq:lagrange}
\end{equation}

\begin{tikzpicture}[scale=2]
  \draw [<->, thick] (0, 3.25) |- (4.25, 0);
  \foreach \x/\y [count=\i] in {2/1, 3/0.5, 4/1.5} {
    \fill (\x, \y) circle [radius=1pt];
    \draw (\x, 0.1) -- (\x, -0.1) node [below] {$x_\i$};
    \draw (0.1, \y) -- (-0.1, \y) node [left] {$y_\i$};
  }
  \draw [densely dotted, gray] (1, 3) circle [radius=1pt];
  \draw [gray] (1, 0.1) -- ++(down:0.2) node [below] {$x_0$};
  \draw [gray] (0.1, 3) -- ++(left:0.2) node [left] {$y_0$};
  \draw [gray, densely dotted] (1, 3) plot [domain=1:4] ({\x}, {0.75*\x^2 - 4.25*\x + 6.5});
  
  \draw [blue] (2.25, 47/64) circle [radius=1.5pt];
  \draw [blue] (2.25, 0.1) -- ++(down:0.2) node [below] {$x_i$};
  \draw [blue] (0.1, 47/64) -- ++(left:0.2) node [left] {$y_i$};
  
  \draw [gray, |<->|] (2, 0.1) -- (2.25, 0.1) node [above, midway] {$h$};
\end{tikzpicture}

One way to interpolate is to pad with the extrapolated values $y_0 = q(x_0) = 3 \, y_1 - 3 \, y_2 + y_3$ and then use the cubic convolutional algorithm as usual to interpolate $y_i$. For a 2D array \lstinline|A|, a $5 \times 5$ matrix in the following example, this could be done in two steps:\\
\begin{minipage}{0.52\linewidth}
  \begin{lstlisting}[gobble=4]
    % Add rows at top and bottom, yields 7x5
    B = [3*A(1, :) - 3*A(2, :) + A(3, :);
        A;
        3*A(end-1, :) - 3*A(end-2, :) + A(end-3, :)];
        
    % Add columns left and right, yields 7x7
    C = [3*B(:, 1) - 3*B(:, 2) + B(:, 3), ...
        B, ...
        3*B(:, end-1) - 3*B(:, end-2) + B(:, end-3)];

  \end{lstlisting}
  This can be combined into one step to save some memory. Even more memory could be saved by interpolating the borders separately with $N \times 4$ and $4 \times M$ padded matrices at the cost of code complexity.
  
  \parindent\keptparindent\ignorespaces
  When the interpolation coordinate match exactly the coordinate of a row or column, it only has to be copied. The kernel function can also handle this, since $h = 0 \Rightarrow W(h) = 1$ for these cases. However, in the special case of copying the last row or column, there are missing two neighbors again. This can be caught and handled separately or solved by padding an additional zero row and column.
\end{minipage}%
\hfill%
\begin{minipage}{0.45\linewidth}
  \centering
  \begin{tikzpicture}[yscale=-1, x=12mm, y=12mm, font=\ttfamily]
    \foreach \y in {1, ..., 5} {
      \foreach \x in {1, ..., 5} {
        \fill (\x, \y) circle [radius=1.5pt];
      }
    }
    \node at (3, 3.5) {A};
    \node [red] at (3, 0.5) {B};
    \node [blue] at (0.5, 3) {C};
    \foreach \x in {1, ..., 5} {
      \fill [red] (\x, 0) circle [radius=1.5pt];
      \fill [red] (\x, 6) circle [radius=1.5pt];
    }
    \foreach \y in {0, ..., 6} {
      \fill [blue] (0, \y) circle [radius=1.5pt];
      \fill [blue] (6, \y) circle [radius=1.5pt];
    }
  \end{tikzpicture}
\end{minipage}
\bigskip

The second way is to derive a quadratic convolutional kernel from \eqref{eq:lagrange} and apply it directly without padding. This saves memory, but increases code complexity a lot. It holds
\begin{align*}
  y_i =
  q(x_i) &= \frac 1 2 \, (h - 1) \, (h - 2) \, y_1
          - h \, (h - 2) \, y_2
          + \frac 1 2 \, h  \, (h - 1) \, y_3\\
         &= {\underbrace{\left( \frac 1 2 \, h^2 - 3 \, h + 2 \right)}_{(I)}} y_1
          + {\underbrace{\left( -h^2 + 2 \, h \right)}_{(II)}} y_2
          + {\underbrace{\left( \frac 1 2 \, h^2 - \frac 1 2 \, h \right)}_{(III)}} y_3.
\end{align*}
We want to use $(I)$ with $-h$ on $[-1, 0)$, $(II)$ with $1-h$ on $[0, 1)$ and $(III)$ with $2-h$ on $[1, 2)$. Therefore the quadratic kernel function is\\
\begin{minipage}{0.38\linewidth}
  \begin{tikzpicture}[x=9mm, y=9mm]
    \draw [->, thick] (-3, 0) -- (3, 0) node [below] {$h$};
    \draw [->, thick] (0, -0.15) -- (0, 1.25) node [right, blue] {$V_l(h)$};
    \draw [gray] (0.1, 1) -- (-0.1, 1) node [left] {$1$};
    \foreach \x in {-1, ..., 2}
      \draw [gray] (\x, 0.1) -- (\x, -0.1) node [below] {$\x$};
    \draw [blue] (-3, 0) -- (-1, 0)
                         plot [domain=-1:0] ({\x}, {0.5 * \x*\x + 1.5 * \x + 1})
                         plot [domain=0:1] ({\x}, {-\x*\x + 1})
                         plot [domain=1:2] ({\x}, {0.5 * \x*\x - 1.5 * \x + 1})
                         -- (3, 0);
  \end{tikzpicture}
\end{minipage}%
\begin{minipage}{0.62\linewidth}
  \begin{equation}
    V_l(h) =
    \begin{cases}
      \frac 1 2 \, h^2 + \frac 3 2 \, h + 1 & \text{if } -1 \le h < 0\\
      -h^2 + 1 & \text{if } 0 \le h < 1\\
      \frac 1 2 \, h^2 - \frac 3 2 \, h + 1 & \text{if } 1 \le h < 2\\
      0 & \text{otherwise}
    \end{cases}
  \end{equation}
\end{minipage}%

\noindent
For the right border it holds $h = x_i - x_2$, since $x_i$ is between $x_2$ and $x_3$. Hence, we need a slightly different kernel function.
\begin{align*}
  y_i =
  q(x_i) &= {\underbrace{\left( \frac 1 2 \, h^2 + \frac 1 2 \, h \right)}_{(I)}} y_1
          + {\underbrace{\left( -h^2 + 1 \right)}_{(II)}} y_2
          + {\underbrace{\left( \frac 1 2 \, h^2 - \frac 1 2 \, h \right)}_{(III)}} y_3.
\end{align*}
We want to use $(I)$ with $-1-h$ on $[-2, -1)$, $(II)$ with $0-h$ on $[-1, 0)$ and $(III)$ with $1-h$ on $[0, 1)$. Therefore the quadratic kernel function for the right border is\\
\begin{minipage}{0.38\linewidth}
  \begin{tikzpicture}[x=9mm, y=9mm]
    \draw [->, thick] (-3, 0) -- (3, 0) node [below] {$h$};
    \draw [->, thick] (0, -0.15) -- (0, 1.25) node [right, blue] {$V_r(h)$};
    \draw [gray] (0.1, 1) -- (-0.1, 1) node [left] {$1$};
    \foreach \x in {-2, ..., 1}
      \draw [gray] (\x, 0.1) -- (\x, -0.1) node [below] {$\x$};
    \draw [blue] (-3, 0) -- (-2, 0)
                         plot [domain=-2:-1] ({\x}, {0.5 * \x*\x + 1.5 * \x + 1})
                         plot [domain=-1:0] ({\x}, {- \x*\x + 1})
                         plot [domain=0:1] ({\x}, {0.5 * \x*\x - 1.5 * \x + 1})
                         -- (3, 0);
  \end{tikzpicture}
\end{minipage}%
\begin{minipage}{0.62\linewidth}
  \begin{equation}
    V_r(h) =
    \begin{cases}
      \frac 1 2 \, h^2 + \frac 3 2 \, h + 1 & \text{if } -2 \le h < -1\\
      -h^2 + 1 & \text{if } -1 \le h < 0\\
      \frac 1 2 \, h^2 - \frac 3 2 \, h + 1 & \text{if } 0 \le h < 1\\
      0 & \text{otherwise}
    \end{cases}
  \end{equation}
\end{minipage}%

\end{preview}
\end{document}