Бикуби́ческая интерполя́ция — в вычислительной математике расширение кубической интерполяции на случай функции двух переменных, значения которой заданы на двумерной регулярной сетке. Поверхность, полученная в результате бикубической интерполяции, является гладкой функцией на границах соседних квадратов, в отличие от поверхностей, полученных в результате билинейной интерполяции или интерполяции методом ближайшего соседа .
Результат бикубической интерполяции функции заданной в ячейке 4х4
[
0
,
3
]
×
[
0
,
3
]
{\displaystyle [0,3]\times [0,3]}
Сравнение разных методов одномерной и двумерной интерполяций Бикубическая интерполяция часто используется в обработке изображений , давая более качественную картинку по сравнению с билинейной интерполяцией. Также бикубическая интерполяция применяется в алгоритмах управления станков с ЧПУ для учёта неровностей плоскостей, например, при фрезеровке печатных плат.
Принцип метода
править
В случае бикубической интерполяции значение функции в искомой точке вычисляется через её значения в 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]}
![{\displaystyle [0,3]\times [0,3]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/22b08d76f2699e7a750a0af9c1cba8a6e09995c3)
![{\displaystyle [0,1]\times [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/92f35a051af39d8299688d7c4a63e39ee5f95c8b)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/006af2a6715fb86137b706f259e09b884d9527dc)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/23b8e31888ed2c3635b625d5ec47a51de2777455)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3deb232059319ea61bedc3f46099cf43f8d628a8)
![{\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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0f579294aea235a7cf7b1bcabfa56c04df1e2d5d)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/df6a2a1081f09cd78794339e043becbc5b202af1)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1083746c4109786271d1995e32d1c812e394263c)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6b896e8fba67cb2b7220533cc1f708ede37bfe1e)
![{\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}=}](https://wikimedia.org/api/rest_v1/media/math/render/svg/253b350bb03b17b9dc3c6c9a45156ccab551b012)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c34c08ec1ddcfa2c27d518693a84ac8842867325)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a6bf843bfe8106881ab8ecd2bfe26d14b813ca67)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7577f90259b8d3345fd02a07765b1a708afd1589)