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

Расширение кубической сплайн-интерполяции

Сравнение бикубической интерполяции с некоторыми 1- и 2-мерными интерполяциями.
Черные и красные / желтые / зеленые / синие точки соответствуют интерполированной точке и соседним образцам соответственно.
Их высоты над землей соответствуют их значениям.

В математике бикубическая интерполяция является расширением кубической сплайн-интерполяции (метода применения кубической интерполяции к набору данных) для интерполяции точек данных на двумерной регулярной сетке . Интерполированная поверхность (имеется в виду форма ядра, а не изображение) более гладкая , чем соответствующие поверхности, полученные билинейной интерполяцией или интерполяцией по ближайшим соседям . Бикубическая интерполяция может быть выполнена с использованием либо полиномов Лагранжа , либо кубических сплайнов , либо алгоритма кубической свертки.

При обработке изображений бикубическая интерполяция часто выбирается вместо билинейной или интерполяции ближайшего соседа при повторной выборке изображения , когда скорость не имеет значения. В отличие от билинейной интерполяции, которая учитывает только 4 пикселя (2×2), бикубическая интерполяция учитывает 16 пикселей (4×4). Изображения, повторно выбранные с помощью бикубической интерполяции, могут иметь различные артефакты интерполяции в зависимости от выбранных значений b и c.

Вычисление

Бикубическая интерполяция на квадрате, состоящем из 25 единичных квадратов, соединенных вместе. Бикубическая интерполяция согласно реализации Matplotlib . Цвет указывает значение функции. Черные точки — это местоположения заданных интерполируемых данных. Обратите внимание, что цветовые образцы не являются радиально-симметричными. ${\displaystyle [0,4]\times [0,4]}$

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

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

Предположим, что значения функции и производные , и известны в четырех углах , , , и единичного квадрата. Тогда интерполированную поверхность можно записать как ${\displaystyle f}$ ${\displaystyle f_{x}}$ ${\displaystyle f_{y}}$ ${\displaystyle f_{xy}}$ ${\displaystyle (0,0)}$ ${\displaystyle (1,0)}$ ${\displaystyle (0,1)}$ ${\displaystyle (1,1)}$ $${\displaystyle p(x,y)=\sum \limits _{i=0}^{3}\sum _{j=0}^{3}a_{ij}x^{i}y^{j}.}$$

Задача интерполяции состоит в определении 16 коэффициентов . Сопоставление со значениями функции дает четыре уравнения: ${\displaystyle a_{ij}}$ ${\displaystyle p(x,y)}$

1. ${\displaystyle f(0,0)=p(0,0)=a_{00},}$
2. ${\displaystyle f(1,0)=p(1,0)=a_{00}+a_{10}+a_{20}+a_{30},}$
3. ${\displaystyle f(0,1)=p(0,1)=a_{00}+a_{01}+a_{02}+a_{03},}$
4. ${\displaystyle f(1,1)=p(1,1)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=0}^{3}a_{ij}.}$

Аналогично, восемь уравнений для производных по направлениям и : ${\displaystyle x}$ ${\displaystyle y}$

1. ${\displaystyle f_{x}(0,0)=p_{x}(0,0)=a_{10},}$
2. ${\displaystyle f_{x}(1,0)=p_{x}(1,0)=a_{10}+2a_{20}+3a_{30},}$
3. ${\displaystyle f_{x}(0,1)=p_{x}(0,1)=a_{10}+a_{11}+a_{12}+a_{13},}$
4. ${\displaystyle f_{x}(1,1)=p_{x}(1,1)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}a_{ij}i,}$
5. ${\displaystyle f_{y}(0,0)=p_{y}(0,0)=a_{01},}$
6. ${\displaystyle f_{y}(1,0)=p_{y}(1,0)=a_{01}+a_{11}+a_{21}+a_{31},}$
7. ${\displaystyle f_{y}(0,1)=p_{y}(0,1)=a_{01}+2a_{02}+3a_{03},}$
8. ${\displaystyle f_{y}(1,1)=p_{y}(1,1)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=1}^{3}a_{ij}j.}$

И четыре уравнения для смешанной частной производной : ${\displaystyle xy}$

1. ${\displaystyle f_{xy}(0,0)=p_{xy}(0,0)=a_{11},}$
2. ${\displaystyle f_{xy}(1,0)=p_{xy}(1,0)=a_{11}+2a_{21}+3a_{31},}$
3. ${\displaystyle f_{xy}(0,1)=p_{xy}(0,1)=a_{11}+2a_{12}+3a_{13},}$
4. ${\displaystyle f_{xy}(1,1)=p_{xy}(1,1)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=1}^{3}a_{ij}ij.}$

В приведенных выше выражениях использованы следующие идентичности: $${\displaystyle p_{x}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}a_{ij}ix^{i-1}y^{j},}$$ $${\displaystyle p_{y}(x,y)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=1}^{3}a_{ij}x^{i}jy^{j-1},}$$ $${\displaystyle p_{xy}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=1}^{3}a_{ij}ix^{i-1}jy^{j-1}.}$$

Эта процедура дает поверхность на единичном квадрате , которая является непрерывной и имеет непрерывные производные. Бикубическая интерполяция на произвольной размерной регулярной сетке может быть затем выполнена путем склеивания таких бикубических поверхностей, гарантируя, что производные совпадают на границах. ${\displaystyle p(x,y)}$ ${\displaystyle [0,1]\times [0,1]}$

Группируем неизвестные параметры в вектор и позволяем приведенной выше системе уравнений быть переформулированной в матрицу для линейного уравнения . ${\displaystyle a_{ij}}$ $${\displaystyle \alpha =\left[{\begin{smallmatrix}a_{00}&a_{10}&a_{20}&a_{30}&a_{01}&a_{11}&a_{21}&a_{31}&a_{02}&a_{12}&a_{22}&a_{32}&a_{03}&a_{13}&a_{23}&a_{33}\end{smallmatrix}}\right]^{T}}$$ $${\displaystyle x=\left[{\begin{smallmatrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&f_{x}(0,0)&f_{x}(1,0)&f_{x}(0,1)&f_{x}(1,1)&f_{y}(0,0)&f_{y}(1,0)&f_{y}(0,1)&f_{y}(1,1)&f_{xy}(0,0)&f_{xy}(1,0)&f_{xy}(0,1)&f_{xy}(1,1)\end{smallmatrix}}\right]^{T},}$$ ${\displaystyle A\alpha =x}$

Обращение матрицы дает более полезное линейное уравнение , где которое позволяет быстро и легко производить вычисления. ${\displaystyle A^{-1}x=\alpha }$ $${\displaystyle A^{-1}=\left[{\begin{smallmatrix}{\begin{array}{rrrrrrrrrrrrrrrr}1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\-3&3&0&0&-2&-1&0&0&0&0&0&0&0&0&0&0\\2&-2&0&0&1&1&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&1&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0&0&0&1&0&0&0\\0&0&0&0&0&0&0&0&-3&3&0&0&-2&-1&0&0\\0&0&0&0&0&0&0&0&2&-2&0&0&1&1&0&0\\-3&0&3&0&0&0&0&0&-2&0&-1&0&0&0&0&0\\0&0&0&0&-3&0&3&0&0&0&0&0&-2&0&-1&0\\9&-9&-9&9&6&3&-6&-3&6&-6&3&-3&4&2&2&1\\-6&6&6&-6&-3&-3&3&3&-4&4&-2&2&-2&-2&-1&-1\\2&0&-2&0&0&0&0&0&1&0&1&0&0&0&0&0\\0&0&0&0&2&0&-2&0&0&0&0&0&1&0&1&0\\-6&6&6&-6&-4&-2&4&2&-3&3&-3&3&-2&-1&-2&-1\\4&-4&-4&4&2&2&-2&-2&2&-2&2&-2&1&1&1&1\end{array}}\end{smallmatrix}}\right],}$$ ${\displaystyle \alpha }$

Может быть и другая краткая матричная форма для 16 коэффициентов: или где $${\displaystyle {\begin{bmatrix}f(0,0)&f(0,1)&f_{y}(0,0)&f_{y}(0,1)\\f(1,0)&f(1,1)&f_{y}(1,0)&f_{y}(1,1)\\f_{x}(0,0)&f_{x}(0,1)&f_{xy}(0,0)&f_{xy}(0,1)\\f_{x}(1,0)&f_{x}(1,1)&f_{xy}(1,0)&f_{xy}(1,1)\end{bmatrix}}={\begin{bmatrix}1&0&0&0\\1&1&1&1\\0&1&0&0\\0&1&2&3\end{bmatrix}}{\begin{bmatrix}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{bmatrix}}{\begin{bmatrix}1&1&0&0\\0&1&1&1\\0&1&0&2\\0&1&0&3\end{bmatrix}},}$$ $${\displaystyle {\begin{bmatrix}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{bmatrix}}={\begin{bmatrix}1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1\end{bmatrix}}{\begin{bmatrix}f(0,0)&f(0,1)&f_{y}(0,0)&f_{y}(0,1)\\f(1,0)&f(1,1)&f_{y}(1,0)&f_{y}(1,1)\\f_{x}(0,0)&f_{x}(0,1)&f_{xy}(0,0)&f_{xy}(0,1)\\f_{x}(1,0)&f_{x}(1,1)&f_{xy}(1,0)&f_{xy}(1,1)\end{bmatrix}}{\begin{bmatrix}1&0&-3&2\\0&0&3&-2\\0&1&-2&1\\0&0&-1&1\end{bmatrix}},}$$ $${\displaystyle p(x,y)={\begin{bmatrix}1&x&x^{2}&x^{3}\end{bmatrix}}{\begin{bmatrix}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{bmatrix}}{\begin{bmatrix}1\\y\\y^{2}\\y^{3}\end{bmatrix}}.}$$

Расширение до прямолинейных сеток

Часто приложения требуют бикубической интерполяции с использованием данных на прямоугольной сетке, а не единичного квадрата. В этом случае тождества для и становятся где — интервал ячейки, содержащей точку , и аналогичными для . В этом случае наиболее практичным подходом к вычислению коэффициентов является то, чтобы позволить затем решить с помощью , как и прежде. Затем нормализованные интерполирующие переменные вычисляются как где и — и координаты точек сетки, окружающих точку . Затем интерполирующая поверхность становится ${\displaystyle p_{x},p_{y},}$ ${\displaystyle p_{xy}}$ $${\displaystyle p_{x}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}{\frac {a_{ij}ix^{i-1}y^{j}}{\Delta x}},}$$ $${\displaystyle p_{y}(x,y)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=1}^{3}{\frac {a_{ij}x^{i}jy^{j-1}}{\Delta y}},}$$ $${\displaystyle p_{xy}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=1}^{3}{\frac {a_{ij}ix^{i-1}jy^{j-1}}{\Delta x\Delta y}},}$$ ${\displaystyle \Delta x}$ ${\displaystyle x}$ ${\displaystyle (x,y)}$ ${\displaystyle \Delta y}$ ${\displaystyle \alpha }$ $${\displaystyle x=\left[{\begin{smallmatrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&\Delta xf_{x}(0,0)&\Delta xf_{x}(1,0)&\Delta xf_{x}(0,1)&\Delta xf_{x}(1,1)&\Delta yf_{y}(0,0)&\Delta yf_{y}(1,0)&\Delta yf_{y}(0,1)&\Delta yf_{y}(1,1)&\Delta x\Delta yf_{xy}(0,0)&\Delta x\Delta yf_{xy}(1,0)&\Delta x\Delta yf_{xy}(0,1)&\Delta x\Delta yf_{xy}(1,1)\end{smallmatrix}}\right]^{T},}$$ ${\displaystyle \alpha =A^{-1}x}$ ${\displaystyle A}$ $${\displaystyle {\begin{aligned}{\overline {x}}&={\frac {x-x_{0}}{x_{1}-x_{0}}},\\{\overline {y}}&={\frac {y-y_{0}}{y_{1}-y_{0}}}\end{aligned}}}$$ ${\displaystyle x_{0},x_{1},y_{0},}$ ${\displaystyle y_{1}}$ ${\displaystyle x}$ ${\displaystyle y}$ ${\displaystyle (x,y)}$ $${\displaystyle p(x,y)=\sum \limits _{i=0}^{3}\sum _{j=0}^{3}a_{ij}{\overline {x}}^{i}{\overline {y}}^{j}.}$$

Нахождение производных от значений функции

Если производные неизвестны, их обычно аппроксимируют по значениям функции в точках, соседних с углами единичного квадрата, например, с помощью конечных разностей .

Чтобы найти одну из отдельных производных или , используя этот метод, найдите наклон между двумя окружающими точками на соответствующей оси. Например, чтобы вычислить для одной из точек, найдите для точек слева и справа от целевой точки и вычислите их наклон, и аналогично для . ${\displaystyle f_{x}}$ ${\displaystyle f_{y}}$ ${\displaystyle f_{x}}$ ${\displaystyle f(x,y)}$ ${\displaystyle f_{y}}$

Чтобы найти перекрестную производную , возьмите производную по обеим осям, по одной за раз. Например, можно сначала использовать процедуру для нахождения производных точек выше и ниже целевой точки, затем использовать процедуру на этих значениях (а не, как обычно, значениях для этих точек), чтобы получить значение для целевой точки. (Или можно сделать это в обратном направлении, сначала вычислив , а затем из них. Оба дают эквивалентные результаты.) ${\displaystyle f_{xy}}$ ${\displaystyle f_{x}}$ ${\displaystyle x}$ ${\displaystyle f_{y}}$ ${\displaystyle f}$ ${\displaystyle f_{xy}(x,y)}$ ${\displaystyle f_{y}}$ ${\displaystyle f_{x}}$

На краях набора данных, когда отсутствуют некоторые окружающие точки, недостающие точки могут быть аппроксимированы несколькими методами. Простой и распространенный метод заключается в том, чтобы предположить, что наклон от существующей точки до целевой продолжается без дальнейших изменений, и использовать это для вычисления гипотетического значения для недостающей точки.

Алгоритм бикубической свертки

Бикубическая сплайн-интерполяция требует решения линейной системы, описанной выше, для каждой ячейки сетки. Интерполятор с аналогичными свойствами может быть получен путем применения свертки со следующим ядром в обоих измерениях: где обычно устанавливается равным −0,5 или −0,75. Обратите внимание, что и для всех ненулевых целых чисел . $${\displaystyle W(x)={\begin{cases}(a+2)|x|^{3}-(a+3)|x|^{2}+1&{\text{for }}|x|\leq 1,\\a|x|^{3}-5a|x|^{2}+8a|x|-4a&{\text{for }}1<|x|<2,\\0&{\text{otherwise}},\end{cases}}}$$ ${\displaystyle a}$ ${\displaystyle W(0)=1}$ ${\displaystyle W(n)=0}$ ${\displaystyle n}$

Этот подход был предложен Кейсом, который показал, что он обеспечивает сходимость третьего порядка относительно интервала выборки исходной функции. ^[1] ${\displaystyle a=-0.5}$

Если мы используем матричную запись для общего случая , мы можем выразить уравнение более удобным образом: для между 0 и 1 для одного измерения. Обратите внимание, что для интерполяции одномерной кубической свертки требуются 4 точки выборки. Для каждого запроса два образца расположены слева и два образца справа. Эти точки индексируются от −1 до 2 в этом тексте. Расстояние от точки с индексом 0 до точки запроса обозначается здесь . ${\displaystyle a=-0.5}$ $${\displaystyle p(t)={\tfrac {1}{2}}{\begin{bmatrix}1&t&t^{2}&t^{3}\end{bmatrix}}{\begin{bmatrix}0&2&0&0\\-1&0&1&0\\2&-5&4&-1\\-1&3&-3&1\end{bmatrix}}{\begin{bmatrix}f_{-1}\\f_{0}\\f_{1}\\f_{2}\end{bmatrix}}}$$ ${\displaystyle t}$ ${\displaystyle t}$

Для двух измерений, впервые примененных один раз в и еще раз в : ${\displaystyle x}$ ${\displaystyle y}$ $${\displaystyle {\begin{aligned}b_{-1}&=p(t_{x},f_{(-1,-1)},f_{(0,-1)},f_{(1,-1)},f_{(2,-1)}),\\[1ex]b_{0}&=p(t_{x},f_{(-1,0)},f_{(0,0)},f_{(1,0)},f_{(2,0)}),\\[1ex]b_{1}&=p(t_{x},f_{(-1,1)},f_{(0,1)},f_{(1,1)},f_{(2,1)}),\\[1ex]b_{2}&=p(t_{x},f_{(-1,2)},f_{(0,2)},f_{(1,2)},f_{(2,2)}),\end{aligned}}}$$ $${\displaystyle p(x,y)=p(t_{y},b_{-1},b_{0},b_{1},b_{2}).}$$

Использование в компьютерной графике

Нижняя половина этого рисунка — это увеличение верхней половины, показывающее, как создается кажущаяся резкость левой линии. Бикубическая интерполяция вызывает перерегулирование, что увеличивает резкость .

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

Однако из-за отрицательных лепестков на ядре это вызывает перерегулирование (ореол). Это может вызвать обрезание и является артефактом (см. также артефакты звона ), но это увеличивает остроту (кажущуюся резкость) и может быть желательным.

Смотрите также

• Математический портал

• Пространственное сглаживание
• Поверхность Безье
• Билинейная интерполяция
• Кубический сплайн Эрмита , одномерный аналог бикубического сплайна
• Повторная выборка Ланцоша
• Интерполяция естественного соседа
• Фильтр синк
• Сплайн-интерполяция
• Трикубическая интерполяция
• Направленная кубическая сверточная интерполяция

Ссылки

1. R. Keys (1981). «Интерполяция кубической свертки для цифровой обработки изображений». Труды IEEE по акустике, речи и обработке сигналов . 29 (6): 1153–1160. Bibcode :1981ITASS..29.1153K. CiteSeerX  10.1.1.320.776 . doi :10.1109/TASSP.1981.1163711.

Внешние ссылки

• Применение интерполяции к выборкам высот
• Теория интерполяции
• Объяснение и реализация (би)кубической интерполяции на Java/C++
• Функция Excel Worksheet для бикубической интерполяции Лагранжа