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