Результат бикубической интерполяции функции заданной в ячейке 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]}
.