Результат бикубической интерполяции функции заданной в ячейке 4х4 [ 0 , 3 ] × [ 0 , 3 ] {\displaystyle [0,3]\times [0,3]} . Данную ячейку можно рассматривать как квадрат, состоящий из 9 единичных квадратов. Черными точками обозначены известные значения функции до интерполяции. Условным цветом обозначены интерполированные значения в каждой точке полученной функции. Сравнение разных методов одномерной и двумерной интерполяций Бикуби́ческая интерполя́ция — в вычислительной математике расширение кубической интерполяции на случай функции двух переменных, значения которой заданы на двумерной регулярной сетке. Поверхность, полученная в результате бикубической интерполяции, является гладкой функцией на границах соседних квадратов, в отличие от поверхностей, полученных в результате билинейной интерполяции или интерполяции методом ближайшего соседа .
Бикубическая интерполяция часто используется в обработке изображений , давая более качественную картинку по сравнению с билинейной интерполяцией. Также бикубическая интерполяция применяется в алгоритмах управления станков с ЧПУ для учёта неровностей плоскостей, например, при фрезеровке печатных плат.
В случае бикубической интерполяции значение функции в искомой точке вычисляется через её значения в 16 соседних точках, расположенных в вершинах квадратов плоскости x , y {\displaystyle x,y} .
При использовании приведённых ниже формул для программной реализации бикубической интерполяции следует помнить, что значения x {\displaystyle x} и y {\displaystyle y} являются относительными, а не абсолютными. Например, для точки с координатами ( 100.3 , 100.8 ) {\displaystyle (100.3,100.8)} x = 0.3 , y = 0.8 {\displaystyle 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 ) + {\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)+} + 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_{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 = 1 4 ( x − 1 ) ( x − 2 ) ( x + 1 ) ( y − 1 ) ( y − 2 ) ( y + 1 ) {\displaystyle b_{1}={\frac {1}{4}}(x-1)(x-2)(x+1)(y-1)(y-2)(y+1)} , b 2 = − 1 4 x ( x + 1 ) ( x − 2 ) ( y − 1 ) ( y − 2 ) ( y + 1 ) {\displaystyle b_{2}=-{\frac {1}{4}}x(x+1)(x-2)(y-1)(y-2)(y+1)} , b 3 = − 1 4 y ( x − 1 ) ( x − 2 ) ( x + 1 ) ( y + 1 ) ( y − 2 ) {\displaystyle b_{3}=-{\frac {1}{4}}y(x-1)(x-2)(x+1)(y+1)(y-2)} , b 4 = 1 4 x y ( x + 1 ) ( x − 2 ) ( y + 1 ) ( y − 2 ) {\displaystyle b_{4}={\frac {1}{4}}xy(x+1)(x-2)(y+1)(y-2)} , b 5 = − 1 12 x ( x − 1 ) ( x − 2 ) ( y − 1 ) ( y − 2 ) ( y + 1 ) {\displaystyle b_{5}=-{\frac {1}{12}}x(x-1)(x-2)(y-1)(y-2)(y+1)} , b 6 = − 1 12 y ( x − 1 ) ( x − 2 ) ( x + 1 ) ( y − 1 ) ( y − 2 ) {\displaystyle b_{6}=-{\frac {1}{12}}y(x-1)(x-2)(x+1)(y-1)(y-2)} , b 7 = 1 12 x y ( x − 1 ) ( x − 2 ) ( y + 1 ) ( y − 2 ) {\displaystyle b_{7}={\frac {1}{12}}xy(x-1)(x-2)(y+1)(y-2)} , b 8 = 1 12 x y ( x + 1 ) ( x − 2 ) ( y − 1 ) ( y − 2 ) {\displaystyle b_{8}={\frac {1}{12}}xy(x+1)(x-2)(y-1)(y-2)} , b 9 = 1 12 x ( x − 1 ) ( x + 1 ) ( y − 1 ) ( y − 2 ) ( y + 1 ) {\displaystyle b_{9}={\frac {1}{12}}x(x-1)(x+1)(y-1)(y-2)(y+1)} , b 10 = 1 12 y ( x − 1 ) ( x − 2 ) ( x + 1 ) ( y − 1 ) ( y + 1 ) {\displaystyle b_{10}={\frac {1}{12}}y(x-1)(x-2)(x+1)(y-1)(y+1)} , b 11 = 1 36 x y ( x − 1 ) ( x − 2 ) ( y − 1 ) ( y − 2 ) {\displaystyle b_{11}={\frac {1}{36}}xy(x-1)(x-2)(y-1)(y-2)} , b 12 = − 1 12 x y ( x − 1 ) ( x + 1 ) ( y + 1 ) ( y − 2 ) {\displaystyle b_{12}=-{\frac {1}{12}}xy(x-1)(x+1)(y+1)(y-2)} , b 13 = − 1 12 x y ( x + 1 ) ( x − 2 ) ( y − 1 ) ( y + 1 ) {\displaystyle b_{13}=-{\frac {1}{12}}xy(x+1)(x-2)(y-1)(y+1)} , b 14 = − 1 36 x y ( x − 1 ) ( x + 1 ) ( y − 1 ) ( y − 2 ) {\displaystyle b_{14}=-{\frac {1}{36}}xy(x-1)(x+1)(y-1)(y-2)} , b 15 = − 1 36 x y ( x − 1 ) ( x − 2 ) ( y − 1 ) ( y + 1 ) {\displaystyle b_{15}=-{\frac {1}{36}}xy(x-1)(x-2)(y-1)(y+1)} , b 16 = 1 36 x y ( x − 1 ) ( x + 1 ) ( y − 1 ) ( y + 1 ) {\displaystyle b_{16}={\frac {1}{36}}xy(x-1)(x+1)(y-1)(y+1)} , Подобным образом можно использовать и интерполяции более высокого порядка, вычисляя значения функции по соседним 4 k 2 {\displaystyle 4k^{2}} точкам.
Результат билинейной интерполяции на тех же входных данных. Частные производные не являются непрерывными и терпят разрыв на границах квадратов 4х4. Результат интерполяции методом ближайшего соседа на тех же входных данных. Допустим, что необходимо интерполировать значение функции f ( x , y ) {\displaystyle f(x,y)} в точке P ( x , y ) {\displaystyle P(x,y)} , лежащей внутри квадрата [ 0 , 1 ] × [ 0 , 1 ] {\displaystyle [0,1]\times [0,1]} , и известно значение функции f {\displaystyle f} в шестнадцати соседних точках ( i , j ) , i = − 1 … 2 , j = − 1 … 2 {\displaystyle (i,j),i=-1\ldots 2,j=-1\ldots 2} .
Тогда общий вид функции, задающей интерполированную поверхность, может быть записан следующим образом:
p ( x , y ) = ∑ i = 0 3 ∑ j = 0 3 a i j x i y j {\displaystyle p(x,y)=\sum _{i=0}^{3}\sum _{j=0}^{3}a_{ij}x^{i}y^{j}} .
Для нахождения коэффициентов a i j {\displaystyle a_{ij}} необходимо подставить в вышеприведённое уравнение значения функции в известных шестнадцати точках. Например:
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 + 2 a 10 + 4 a 20 + 8 a 30 {\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}}} .
Полностью в матричном виде:
M α T = γ T {\displaystyle M\alpha ^{T}=\gamma ^{T}} ,
где
α = [ 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 ] {\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]} ,
γ = [ 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 ) ] {\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]} ,
M = [ 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 ] {\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]} .
Решая получившуюся систему линейных алгебраических уравнений , можно найти значения a i j {\displaystyle a_{ij}} в явном виде:
α T = 1 36 [ 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 ] γ T {\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}} .
Единожды найденные коэффициенты a i j {\displaystyle a_{ij}} теперь могут быть использованы для многократного вычисления интерполированного значения функции в произвольных точках квадрата [ 0 , 1 ] × [ 0 , 1 ] {\displaystyle [0,1]\times [0,1]} .
Следует отметить, что такой способ обеспечивает непрерывность самой функции и её второй производной на границах смежных квадратов, но приводит к разрыву первых производных на границах ячеек 4×4. Для обеспечения непрерывности самой функции и её первой производной необходимо подставлять в исходное выражение значения функции и значения первых производных по направлениям x и y в вершинах центральной ячейки, производные рассчитываются через центральные разности. Для подстановки производных выражение должно быть продифференцировано соответствующим образом.
Другая интерпретация метода заключается в том, что для нахождения интерполированного значения можно сначала произвести кубическую интерполяцию в одном направлении, а затем в другом.
Для функции f ( x ) {\displaystyle f(x)} с известными значениями f ( − 1 ) {\displaystyle f(-1)} , f ( 0 ) {\displaystyle f(0)} , f ( 1 ) {\displaystyle f(1)} , f ( 2 ) {\displaystyle f(2)} можно построить кубический сплайн: p ( x ) = ∑ i = 0 3 b i x i {\displaystyle p(x)=\textstyle \sum _{i=0}^{3}b_{i}x^{i}} , или в матричном виде:
p ( x ) = [ 1 x x 2 x 3 ] [ b 0 b 1 b 2 b 3 ] {\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]} ,
где
[ b 0 b 1 b 2 b 3 ] = A [ f ( − 1 ) f ( 0 ) f ( 1 ) f ( 2 ) ] {\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]} ,
A = 1 3 [ 0 6 0 0 − 2 − 3 6 − 1 3 − 6 3 0 − 1 3 − 3 1 ] {\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]} .
Таким образом, для нахождения интерполированного значения p ( x , y ) {\displaystyle p(x,y)} в квадрате [ 0 , 1 ] × [ 0 , 1 ] {\displaystyle [0,1]\times [0,1]} можно сначала рассчитать четыре значения p ( x , − 1 ) {\displaystyle p(x,-1)} , p ( x , 0 ) {\displaystyle p(x,0)} , p ( x , 1 ) {\displaystyle p(x,1)} , p ( x , 2 ) {\displaystyle p(x,2)} для зафиксированного x {\displaystyle x} , затем через полученные четыре точки построить кубический сплайн, и этим завершить вычисление p ( x , y ) {\displaystyle p(x,y)} :
p ( x , y ) = [ 1 y y 2 y 3 ] A ( [ 1 x x 2 x 3 ] A [ 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 ) ] ) T = {\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}=}
= [ 1 y y 2 y 3 ] A [ 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 ) ] A T [ 1 x x 2 x 3 ] {\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]} .
Следует отметить, что такой подход обеспечивает непрерывность самой функции и её вторых производных на границе ячеек, но не обеспечивает непрерывности первой производной. Для обеспечения непрерывности первой производной необходимо подставлять значения функции и её первых производных на границе центральной ячейки. Тогда коэффициенты сплайна будут иметь вид:
[ b 0 b 1 b 2 b 3 ] = B − 1 [ f ( 0 ) f ( 1 ) f ( 1 ) − f ( − 1 ) 2 f ( 2 ) − f ( 0 ) 2 ] {\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]} ,
B = [ 1 0 0 0 1 1 1 1 0 1 0 0 0 1 2 3 ] , B − 1 = [ 1 0 0 0 0 0 1 0 − 3 3 − 2 − 1 2 − 2 1 1 ] {\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]} .