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

Результат бикубической интерполяции функции заданной в ячейке 4х4 $[0,3]\times [0,3]$. Данную ячейку можно рассматривать как квадрат, состоящий из 9 единичных квадратов. Черными точками обозначены известные значения функции до интерполяции. Условным цветом обозначены интерполированные значения в каждой точке полученной функции.

Сравнение разных методов одномерной и двумерной интерполяций

Бикуби́ческая интерполя́ция — в вычислительной математике расширение кубической интерполяции на случай функции двух переменных, значения которой заданы на двумерной регулярной сетке. Поверхность, полученная в результате бикубической интерполяции, является гладкой функцией на границах соседних квадратов, в отличие от поверхностей, полученных в результате билинейной интерполяции или интерполяции методом ближайшего соседа.

Бикубическая интерполяция часто используется в обработке изображений, давая более качественную картинку по сравнению с билинейной интерполяцией. Также бикубическая интерполяция применяется в алгоритмах управления станков с ЧПУ для учёта неровностей плоскостей, например, при фрезеровке печатных плат.

Принцип метода

В случае бикубической интерполяции значение функции в искомой точке вычисляется через её значения в 16 соседних точках, расположенных в вершинах квадратов плоскости $x,y$.

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

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

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

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

Допустим, что необходимо интерполировать значение функции $f(x,y)$ в точке $P(x,y)$, лежащей внутри квадрата $[0,1]\times [0,1]$, и известно значение функции $f$ в шестнадцати соседних точках $(i,j),i=-1\dots 2,j=-1\dots 2$.

Тогда общий вид функции, задающей интерполированную поверхность, может быть записан следующим образом:

$p(x,y)=\sum _{i=0}^3\sum _{j=0}^3{a}_{ij}{x}^i{y}^j$.

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

$\begin{array}{cccccccc}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{array}$.

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

$M{\alpha }^T={\gamma }^T$,

где

$\alpha =\left[\begin{array}{cccccccccccccccc}{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{array}\right]$,

$\gamma =\left[\begin{array}{cccccccccccccccc}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{array}\right]$,

$M=\left[\begin{array}{cccccccccccccccc}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{array}\right]$.

Решая получившуюся систему линейных алгебраических уравнений, можно найти значения $a_{ij}$ в явном виде:

${\alpha }^T=\frac{1}{36}\left[\begin{array}{cccccccccccccccc}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{array}\right]{\gamma }^T$.

Единожды найденные коэффициенты $a_{ij}$ теперь могут быть использованы для многократного вычисления интерполированного значения функции в произвольных точках квадрата $[0,1]\times [0,1]$.

Следует отметить, что такой способ обеспечивает непрерывность самой функции и её второй производной на границах смежных квадратов, но приводит к разрыву первых производных на границах ячеек 4×4. Для обеспечения непрерывности самой функции и её первой производной необходимо подставлять в исходное выражение значения функции и значения первых производных по направлениям x и y в вершинах центральной ячейки, производные рассчитываются через центральные разности. Для подстановки производных выражение должно быть продифференцировано соответствующим образом.

Последовательная кубическая интерполяция

Другая интерпретация метода заключается в том, что для нахождения интерполированного значения можно сначала произвести кубическую интерполяцию в одном направлении, а затем в другом.

Для функции $f(x)$ с известными значениями $f(-1)$, $f(0)$, $f(1)$, $f(2)$ можно построить кубический сплайн: $p(x)={\sum }_{i=0}^3{b}_i{x}^i$, или в матричном виде:

$p(x)=\left[\begin{array}{cccc}1& x& x^2& x^3\end{array}\right]\left[\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\end{array}\right]$,

где

$\left[\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\end{array}\right]=A\left[\begin{array}{c}f(-1)\\ f(0)\\ f(1)\\ f(2)\end{array}\right]$,

$A=\frac{1}{3}\left[\begin{array}{cccc}0& 6& 0& 0\\ -2& -3& 6& -1\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\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{array}{cccc}1& y& y^2& y^3\end{array}\right]A{\left(\left[\begin{array}{cccc}1& x& x^2& x^3\end{array}\right]A\left[\begin{array}{ccccc}& 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{array}\right]\right)}^T=$

$=\left[\begin{array}{cccc}1& y& y^2& y^3\end{array}\right]A\left[\begin{array}{cccc}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{array}\right]A^T\left[\begin{array}{c}1\\ x\\ x^2\\ x^3\end{array}\right]$.

Следует отметить, что такой подход обеспечивает непрерывность самой функции и её вторых производных на границе ячеек, но не обеспечивает непрерывности первой производной. Для обеспечения непрерывности первой производной необходимо подставлять значения функции и её первых производных на границе центральной ячейки. Тогда коэффициенты сплайна будут иметь вид:

$\left[\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\end{array}\right]=B^{-1}\left[\begin{array}{c}f(0)\\ f(1)\\ \frac{f(1)-f(-1)}{2}\\ \frac{f(2)-f(0)}{2}\end{array}\right]$,

$B=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& 1& 1\\ 0& 1& 0& 0\\ 0& 1& 2& 3\end{array}\right],B^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ -3& 3& -2& -1\\ 2& -2& 1& 1\end{array}\right]$.

См. также

• Интерполяция
• Билинейная интерполяция
• Фильтр Ланцоша

Литература

• Шаблон:Статья