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

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

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

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

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

${\displaystyle 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)+}$

${\displaystyle +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)}$,

где

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

${\displaystyle {\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}}}$.

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

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

где

${\displaystyle \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]}$,

${\displaystyle \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]}$,

${\displaystyle 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]}$.

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

${\displaystyle \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}}$.

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

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

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

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

${\displaystyle 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]}$,

где

${\displaystyle \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]}$,

${\displaystyle 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]}$.

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

${\displaystyle 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}=}$

${\displaystyle =\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]}$.

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

${\displaystyle \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]}$,

${\displaystyle 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.