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

Материал из Википедии — свободной энциклопедии
Перейти к: навигация, поиск
Результат бикубической интерполяции функции заданной на сетке [0,3] \times [0,3]. Данную сетку можно рассматривать как состоящую из 9 единичных квадратов. Черными точками обозначены известные значения функции до интерполяции. Цветом обозначены интерполированные значения в каждой точке полученного изображения.

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

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

Бикубическая интерполяция сплайнами[править | править вики-текст]

Допустим, что необходимо интерполировать значение функции f(x,y) в точке P(x,y), лежащей внутри квадрата [0,1]\times[0,1], и известно значение функции f в шестнадцати соседних точках (i,j), i=-1\ldots2, j=-1\ldots2. Тогда общий вид функции, задающей интерполированную поверхность, может быть записан следующим образом: p(x,y) = \sum_{i=0}^3 \sum_{j=0}^3 a_{i j} x^i y^j.

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

\begin{align}
f(-1, 0) &=   a_{0 0}  &-  &a_{1 0}  &+  &a_{2 0}  &-  &a_{3 0}  \\
f(0, 0) &=   a_{0 0}  \\
f(1, 0) &=   a_{0 0}  &+  &a_{1 0}  &+  &a_{2 0}  &+  &a_{3 0}  \\
f(2, 0) &=   a_{0 0}  &+ &2  a_{1 0}  &+ &4  a_{2 0}  &+ &8  a_{3 0}  \\
\end{align}.

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

M \alpha^T = \gamma^T,

где

\alpha=\left[\begin{smallmatrix} a_{0 0}  & a_{0 1}  & a_{0 2}  & a_{0 3}  & a_{1 0}  & a_{1 1}  & a_{1 2}  & a_{1 3}  & a_{2 0}  & a_{2 1}  & a_{2 2}  & a_{2 3}  & a_{3 0}  & a_{3 1}  & a_{3 2}  & a_{3 3} \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_{i j} в явном виде:

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

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

См. также[править | править вики-текст]