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

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

Сравнение разных методов одномерной и двумерной интерполяций

Бикуби́ческая интерполя́ция — в вычислительной математике расширение кубической интерполяции на случай функции двух переменных, значения которой заданы на двумерной регулярной сетке. Поверхность, полученная в результате бикубической интерполяции, является гладкой функцией на границах соседних квадратов, в отличие от поверхностей, полученных в результате билинейной интерполяции или интерполяции методом ближайшего соседа.

Бикубическая интерполяция часто используется в обработке изображений, давая более качественную картинку по сравнению с билинейной интерполяцией. Также бикубическая интерполяция применяется в алгоритмах управления станков с ЧПУ для учёта неровностей плоскостей, например, при фрезеровке печатных плат.

Принцип метода

В случае бикубической интерполяции значение функции в искомой точке вычисляется через её значения в 16 соседних точках, расположенных в вершинах квадратов плоскости $x,y$.

При использовании приведённых ниже формул для программной реализации бикубической интерполяции следует помнить, что значения $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}$ точкам.

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

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

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

Допустим, что необходимо интерполировать значение функции $f(x,y)$ в точке $P(x,y)$, лежащей внутри квадрата $[0,1]\times [0,1]$, и известно значение функции $f$ в шестнадцати соседних точках $(i,j),i=-1\ldots 2,j=-1\ldots 2$.

Тогда общий вид функции, задающей интерполированную поверхность, может быть записан следующим образом:

$p(x,y)=\sum _{i=0}^{3}\sum _{j=0}^{3}a_{ij}x^{i}y^{j}$.

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

${\begin{aligned}f(-1,0)&=a_{00}&-&a_{10}&+&a_{20}&-&a_{30}\\f(0,0)&=a_{00}\\f(1,0)&=a_{00}&+&a_{10}&+&a_{20}&+&a_{30}\\f(2,0)&=a_{00}&+&2a_{10}&+&4a_{20}&+&8a_{30}\\\end{aligned}}$.

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

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

где

$\alpha =\left[{\begin{smallmatrix}a_{00}&a_{01}&a_{02}&a_{03}&a_{10}&a_{11}&a_{12}&a_{13}&a_{20}&a_{21}&a_{22}&a_{23}&a_{30}&a_{31}&a_{32}&a_{33}\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_{ij}$ в явном виде:

$\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_{ij}$ теперь могут быть использованы для многократного вычисления интерполированного значения функции в произвольных точках квадрата $[0,1]\times [0,1]$.

Следует отметить, что такой способ обеспечивает непрерывность самой функции и её второй производной на границах смежных квадратов, но приводит к разрыву первых производных на границах ячеек 4×4. Для обеспечения непрерывности самой функции и её первой производной необходимо подставлять в исходное выражение значения функции и значения первых производных по направлениям 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}{3}}\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]$.

См. также

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

Литература

• R. Keys. Cubic convolution interpolation for digital image processing (англ.) // IEEE Transactions on Acoustics, Speech, and Signal Processing : journal. — 1981. — Vol. 29, no. 6. — P. 1153—1160. — doi:10.1109/TASSP.1981.1163711.