Бикубическая интерполяция

Результат бикубической интерполяции функции заданной на сетке $[0,3] \times [0,3]$. Данную сетку можно рассматривать как состоящую из 9 единичных квадратов. Черными точками обозначены известные значения функции до интерполяции. Цветом обозначены интерполированные значения в каждой точке полученного изображения.

Бикубическая интерполяция — в вычислительной математике расширение кубической интерполяции на случай функции двух переменных, значения которой заданы на двумерной регулярной сетке. Поверхность, полученная в результате бикубической интерполяции является гладкой функцией, в отличие от поверхностей, полученных в результате билинейной интерполяции или интерполяции методом ближайшего соседа. Также бикубическая интерполяция часто используется в обработке изображений, давая более качественное изображение по сравнению с билинейной интерполяцией. В случае бикубической интерполяции значение функции в искомой точке вычисляется через ее значения в 16 соседних точках. При использовании приведённых ниже формул для программной реализации бикубической интерполяции следует помнить, что значения $x$ и $y$ являются относительными, а не абсолютными. Например, для пикселя с координатами $(100.3, 100.8)$ $x = 0.3, y = 0.8$. Для получения относительных значений координат необходимо округлить вещественные координаты вниз, и вычесть полученные числа из вещественных координат.

$f(y, x) = b_{1} f(0, 0) + b_{2} f(0, 1) + b_{3} f(1, 0) + b_{4} f(1, 1) + b_{5} f(0, -1) + b_{6} f(-1, 0) + b_{7} f(1, -1) + b_{8} f(-1, 1) +$

$b_{9} f(0, 2) + b_{10} f(2, 0) + b_{11} f(-1,-1) + b_{12} f(1, 2) + b_{13} f(2, 1) + b_{14} f(-1, 2) + b_{15} f(2,-1) + b_{16} f(2, 2)$, где

$b_{1} = \frac{1}{4}(x - 1)(x - 2)(x + 1)(y - 1)(y - 2)(y + 1)$,

$b_{2} = - \frac{1}{4}x(x + 1)(x - 2)(y - 1)(y - 2)(y + 1)$,

$b_{3} = - \frac{1}{4}y(x - 1)(x - 2)(x + 1)(y + 1)(y - 2)$,

$b_{4} = \frac{1}{4}xy(x + 1)(x - 2)(y + 1)(y - 2)$,

$b_{5} = - \frac{1}{12}x(x - 1)(x - 2)(y - 1)(y - 2)(y + 1)$,

$b_{6} = - \frac{1}{12}y(x - 1)(x - 2)(x + 1)(y - 1)(y - 2)$,

$b_{7} = \frac{1}{12}xy(x - 1)(x - 2)(y + 1)(y - 2)$,

$b_{8} = \frac{1}{12}xy(x + 1)(x - 2)(y - 1)(y - 2)$,

$b_{9} = \frac{1}{12}x(x - 1)(x + 1)(y - 1)(y - 2)(y + 1)$,

$b_{10} = \frac{1}{12}y(x - 1)(x - 2)(x + 1)(y - 1)(y + 1)$,

$b_{11} = \frac{1}{36}xy(x - 1)(x - 2)(y - 1)(y - 2)$,

$b_{12} = - \frac{1}{12}xy(x - 1)(x + 1)(y + 1)(y - 2)$,

$b_{13} = - \frac{1}{12}xy(x + 1)(x - 2)(y - 1)(y + 1)$,

$b_{14} = - \frac{1}{36}xy(x - 1)(x + 1)(y - 1)(y - 2)$,

$b_{15} = - \frac{1}{36}xy(x - 1)(x - 2)(y - 1)(y + 1)$,

$b_{16} = \frac{1}{36}xy(x - 1)(x + 1)(y - 1)(y + 1)$,

Подобным образом можно использовать и другие виды интерполяции, вычисляя значения функции по соседним $4k^2$ точкам.

Результат билинейной интерполяции на тех же входных данных. Частные производные не являются непрерывными и терпят разрыв на границах квадратов.

Результат интерполяции методом ближайшего соседа на тех же входных данных.

Бикубическая интерполяция сплайнами

Допустим, что необходимо интерполировать значение функции $f(x,y)$ в точке $P(x,y)$, лежащей внутри квадрата $[0,1]\times[0,1]$, и известно значение функции $f$ в шестнадцати соседних точках $(i,j), i=-1\ldots2, j=-1\ldots2$. Тогда общий вид функции, задающей интерполированную поверхность, может быть записан следующим образом: $p(x,y) = \sum_{i=0}^3 \sum_{j=0}^3 a_{i j} x^i y^j.$

Для нахождения коэффициентов $a_{i j}$ необходимо подставить в выше приведённое уравнение значения функции в известных шестнадцати точках. Например:

$\begin{align} f(-1, 0) &= a_{0 0} &- &a_{1 0} &+ &a_{2 0} &- &a_{3 0} \\ f(0, 0) &= a_{0 0} \\ f(1, 0) &= a_{0 0} &+ &a_{1 0} &+ &a_{2 0} &+ &a_{3 0} \\ f(2, 0) &= a_{0 0} &+ &2 a_{1 0} &+ &4 a_{2 0} &+ &8 a_{3 0} \\ \end{align}$.

Полностью в матричном виде:

$M \alpha^T = \gamma^T$,

где

$\alpha=\left[\begin{smallmatrix} a_{0 0} & a_{0 1} & a_{0 2} & a_{0 3} & a_{1 0} & a_{1 1} & a_{1 2} & a_{1 3} & a_{2 0} & a_{2 1} & a_{2 2} & a_{2 3} & a_{3 0} & a_{3 1} & a_{3 2} & a_{3 3} \end{smallmatrix}\right]$,

$\gamma=\left[\begin{smallmatrix} f(-1, -1) & f(0, -1) & f(1, -1) & f(2, -1) & f(-1, 0) & f(0, 0) & f(1, 0) & f(2, 0) & f(-1, 1) & f(0, 1) & f(1, 1) & f(2, 1) & f(-1, 2) & f(0, 2) & f(1, 2) & f(2, 2) \end{smallmatrix}\right]$,

$M=\left[\begin{smallmatrix} 1 & -1 & 1 & -1 & -1 & 1 & -1 & 1 & 1 & -1 & 1 & -1 & -1 & 1 & -1 & 1 \\ 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 \\ 1 & -1 & 1 & -1 & 2 & -2 & 2 & -2 & 4 & -4 & 4 & -4 & 8 & -8 & 8 & -8 \\ 1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 4 & 0 & 0 & 0 & 8 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 & 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 \\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 2 & 2 & 2 & 2 & 4 & 4 & 4 & 4 & 8 & 8 & 8 & 8 \\ 1 & 2 & 4 & 8 & -1 & -2 & -4 & -8 & 1 & 2 & 4 & 8 & -1 & -2 & -4 & -8 \\ 1 & 2 & 4 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 2 & 4 & 8 & 1 & 2 & 4 & 8 & 1 & 2 & 4 & 8 & 1 & 2 & 4 & 8 \\ 1 & 2 & 4 & 8 & 2 & 4 & 8 & 16 & 4 & 8 & 16 & 32 & 8 & 16 & 32 & 64 \end{smallmatrix}\right]$.

Решая получившуюся систему линейных алгебраических уравнений, можно найти значения $a_{i j}$ в явном виде:

$\alpha^T = \frac{1}{36}\left[\begin{smallmatrix} 0 & 0 & 0 & 0 & 0 & 36 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -12 & 0 & 0 & 0 & -18 & 0 & 0 & 0 & 36 & 0 & 0 & 0 & -6 & 0 & 0 \\ 0 & 18 & 0 & 0 & 0 & -36 & 0 & 0 & 0 & 18 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -6 & 0 & 0 & 0 & 18 & 0 & 0 & 0 & -18 & 0 & 0 & 0 & 6 & 0 & 0 \\ 0 & 0 & 0 & 0 & -12 & -18 & 36 & -6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 4 & 6 & -12 & 2 & 6 & 9 & -18 & 3 & -12 & -18 & 36 & -6 & 2 & 3 & -6 & 1 \\ -6 & -9 & 18 & -3 & 12 & 18 & -36 & 6 & -6 & -9 & 18 & -3 & 0 & 0 & 0 & 0 \\ 2 & 3 & -6 & 1 & -6 & -9 & 18 & -3 & 6 & 9 & -18 & 3 & -2 & -3 & 6 & -1 \\ 0 & 0 & 0 & 0 & 18 & -36 & 18 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -6 & 12 & -6 & 0 & -9 & 18 & -9 & 0 & 18 & -36 & 18 & 0 & -3 & 6 & -3 & 0 \\ 9 & -18 & 9 & 0 & -18 & 36 & -18 & 0 & 9 & -18 & 9 & 0 & 0 & 0 & 0 & 0 \\ -3 & 6 & -3 & 0 & 9 & -18 & 9 & 0 & -9 & 18 & -9 & 0 & 3 & -6 & 3 & 0 \\ 0 & 0 & 0 & 0 & -6 & 18 & -18 & 6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & -6 & 6 & -2 & 3 & -9 & 9 & -3 & -6 & 18 & -18 & 6 & 1 & -3 & 3 & -1 \\ -3 & 9 & -9 & 3 & 6 & -18 & 18 & -6 & -3 & 9 & -9 & 3 & 0 & 0 & 0 & 0 \\ 1 & -3 & 3 & -1 & -3 & 9 & -9 & 3 & 3 & -9 & 9 & -3 & -1 & 3 & -3 & 1 \\ \end{smallmatrix}\right]\gamma^T$.

Единожды найденные коэффициенты $a_{i j}$ теперь могут быть использованы для многократного вычисления интерполированного значения функции в произвольных точках квадрата $[0,1]\times[0,1]$.

Следует отметить, что такой способ обеспечивает непрерывность самой функции и ее второй производной на границах ячеек, но приводит к разрыву первых производных на границах ячеек. Для обеспечения непрерывности самой функции и ее первой производной необходимо подставлять в исходное выражение значения функции и значения первых производных по направлениям x и y в вершинах центральной ячейки, производные рассчитываются через центральные разности. Для подстановки производных выражение должно быть продифферинцировано соответствующим образом.

Последовательная кубическая интерполяция

Другая интерпретация метода заключается в том, что для нахождения интерполированного значения можно сначала произвести кубическую интерполяцию в одном направлении, а затем в другом.

Для функции $f(x)$ с известными значениями $f(-1)$, $f(0)$, $f(1)$, $f(2)$ можно построить кубический сплайн: $p(x) = \textstyle \sum_{i=0}^3 b_i x^i$, или в матричном виде

$p(x) = \left[\begin{smallmatrix} 1 & x & x^2 & x^3 \end{smallmatrix} \right] \left[\begin{smallmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{smallmatrix}\right]$,

где

$\left[\begin{smallmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{smallmatrix}\right] = A \left[\begin{smallmatrix} f(-1) \\ f(0) \\ f(1) \\ f(2) \end{smallmatrix}\right]$,

$A = \frac{1}{6} \left[\begin{smallmatrix} 0 & 6 & 0 & 0\\ -2 & -3 & 6 & -1\\ 3 & -6 & 3 & 0\\ -1 & 3 & -3 & 1 \end{smallmatrix}\right]$.

Таким образом, для нахождения интерполированного значения $p(x,y)$ в квадрате $[0,1]\times[0,1]$ можно сначала рассчитать четыре значения $p(x,-1)$, $p(x,0)$, $p(x,1)$, $p(x,2)$ для зафиксированного $x$, затем через полученные четыре точки построить кубический сплайн, и этим завершить вычисление $p(x,y)$ :

$p(x, y) = \left[\begin{smallmatrix} 1 & y & y^2 & y^3 \end{smallmatrix}\right] A \left( \left[\begin{smallmatrix} 1 & x & x^2 & x^3 \end{smallmatrix}\right] A \left[\begin{smallmatrix} & f(-1, -1) & f(0, -1) & f(1, -1) & f(2, -1) \\ & f(-1, 0) & f(0, 0) & f(1, 0) & f(2, 0) \\ & f(-1, 1) & f(0, 1) & f(1, 1) & f(2, 1) \\ & f(-1, 2) & f(0, 2) & f(1, 2) & f(2, 2) \\ \end{smallmatrix}\right] \right)^T =$

$= \left[\begin{smallmatrix} 1 & y & y^2 & y^3 \end{smallmatrix}\right] A \left[\begin{smallmatrix} f(-1, -1) & f(-1, 0) & f(-1, 1) & f(-1, 2) \\ f(0, -1) & f(0, 0) & f(0, 1) & f(0, 2) \\ f(1, -1) & f(1, 0) & f(1, 1) & f(1, 2) \\ f(2, -1) & f(2, 0) & f(2, 1) & f(2, 2) \end{smallmatrix}\right] A^T \left[\begin{smallmatrix} 1 \\ x \\ x^2 \\ x^3 \end{smallmatrix}\right]$

Следует отметить, что такой подход обеспечивает непрерывность самой функции и ее вторых производных на границе ячеек, но не обеспечивает непрерывности первой производной. Для обеспечения непрерывности первой производной необходимо подставлять значения функции и ее первых производных на границе центральной ячейки. Тогда коэффициенты сплайна будут иметь вид:

$\left[\begin{smallmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{smallmatrix}\right] = B^{-1} \left[\begin{smallmatrix} f(0) \\ f(1) \\ \frac{f(1)-f(-1)}{2} \\ \frac{f(2)-f(0)}{2} \end{smallmatrix}\right]$,

$B = \left[\begin{smallmatrix} 1 & 0 & 0 & 0\\ 1 & 1 & 1 & 1\\ 0 & 1 & 0 & 0\\ 0 & 1 & 2 & 3 \end{smallmatrix}\right], B^{-1} = \left[\begin{smallmatrix} 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0\\ -3 & 3 & -2 & -1\\ 2 & -2 & 1 & 1 \end{smallmatrix}\right]$.

См. также

• Интерполяция
• Билинейная интерполяция
• Фильтр Ланцоша