Interpolación bicúbica

En matemáticas, la interpolación bicúbica es una extensión de la interpolación cúbica (un método para aplicar la interpolación cúbica a un conjunto de datos dispuestos en un retículo regular bidimensional. La superficie interpolada (es decir, la forma del núcleo, no la imagen) es más suave que las superficies correspondientes obtenidas por interpolación bilineal o por interpolación por el vecino más cercano. La interpolación bicúbica se puede lograr utilizando el algoritmo de interpolación polinómica de Lagrange, spline cúbico o convolución cúbica.

En procesamiento digital de imágenes, la interpolación bicúbica a menudo se elige en lugar de la interpolación bilineal o la del vecino más cercano en procesamiento de imágenes, cuando la velocidad del proceso no es un problema. A diferencia de la interpolación bilineal, que solo tiene en cuenta 4 píxeles (2×2), la interpolación bicúbica considera 16 píxeles (4×4). Las imágenes remuestreadas con interpolación bicúbica pueden presentar artefactos de interpolación diferentes, dependiendo de los valores de los parámetros b y c elegidos.

Supóngase que los valores de la función ${\displaystyle f}$ y las derivadas ${\displaystyle f_{x}}$, ${\displaystyle f_{y}}$ y ${\displaystyle f_{xy}}$ se conocen en las cuatro esquinas ${\displaystyle (0,0)}$, ${\displaystyle (1,0)}$, ${\displaystyle (0,1)}$ y ${\displaystyle (1,1)}$ del cuadrado unitario. La superficie interpolada se puede expresar como

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

El problema de interpolación consiste en determinar los 16 coeficientes ${\displaystyle a_{ij}}$. Al hacer coincidir ${\displaystyle p(x,y)}$ con los valores de la función se obtienen cuatro ecuaciones:

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}.}$

Asimismo, se generan ocho ecuaciones para las derivadas en las direcciones ${\displaystyle x}$ e ${\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.}$

Y otras cuatro ecuaciones para la derivada parcial en ${\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.}$

Las expresiones anteriores han utilizado las siguientes identidades:

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

Este procedimiento genera una superficie ${\displaystyle p(x,y)}$ sobre el cuadrado unidad ${\displaystyle [0,1]\times [0,1]}$ que es continua y tiene derivadas continuas. La interpolación bicúbica en un retículo regular de tamaño arbitrario se puede lograr uniendo dichas superficies bicúbicas, asegurando que las derivadas coincidan en los límites.

Agrupando los parámetros desconocidos ${\displaystyle a_{ij}}$ en un vector

${\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}}$

y dejando que

${\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},}$

El sistema de ecuaciones anterior se puede reformular en una matriz para la ecuación lineal ${\displaystyle A\alpha =x}$.

Al invertir la matriz se obtiene la ecuación lineal más manejable ${\displaystyle A^{-1}x=\alpha }$, donde

${\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],}$

lo que permite calcular ${\displaystyle \alpha }$ de forma rápida y sencilla.

Puede haber otra forma matricial concisa para 16 coeficientes:

${\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}},}$

o

${\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}},}$

dónde

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

A menudo, las aplicaciones requieren interpolación bicúbica utilizando datos en una cuadrícula rectilínea, en lugar de la unidad cuadrada. En este caso, las identidades de ${\displaystyle p_{x},p_{y},}$ y ${\displaystyle p_{xy}}$ pasan a ser

${\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}},}$

donde ${\displaystyle \Delta x}$ es el espaciado ${\displaystyle x}$ de la celda que contiene el punto ${\displaystyle (x,y)}$ y similar para ${\displaystyle \Delta y}$. En este caso, el enfoque más práctico para calcular los coeficientes ${\displaystyle \alpha }$ es dejar

${\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},}$

y luego resolver ${\displaystyle \alpha =A^{-1}x}$ con ${\displaystyle A}$ como antes. A continuación, las variables de interpolación normalizadas se calculan como

${\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}}}$

donde ${\displaystyle x_{0},x_{1},y_{0},}$ y ${\displaystyle y_{1}}$ son las coordenadas ${\displaystyle x}$ y ${\displaystyle y}$ de los puntos de la cuadrícula que rodean el punto ${\displaystyle (x,y)}$. Entonces, la superficie de interpolación se convierte en

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

Si se desconocen las derivadas, normalmente se aproximan a partir de los valores de la función en puntos vecinos a las esquinas del cuadrado unitario, utilizando por ejemplo diferencias finitas.

Para encontrar cualquiera de las derivadas simples, ${\displaystyle f_{x}}$ o ${\displaystyle f_{y}}$, usando ese método, se calcula la pendiente entre los dos puntos "circundantes" en el eje apropiado. Por ejemplo, para calcular ${\displaystyle f_{x}}$ para uno de los puntos, se debe buscar ${\displaystyle f(x,y)}$ para los puntos a la izquierda y a la derecha del punto objetivo y calcular su pendiente, y proceder de manera similar para ${\displaystyle f_{y}}$.

Para encontrar la derivada cruzada ${\displaystyle f_{xy}}$, se debe tomar la derivada en ambos ejes, uno cada vez. Por ejemplo, primero se puede usar el procedimiento ${\displaystyle f_{x}}$ para encontrar las derivadas ${\displaystyle x}$ de los puntos por encima y por debajo del punto objetivo, luego usar el procedimiento ${\displaystyle f_{y}}$ en esos valores (en lugar de, como es habitual, los valores de ${\displaystyle f}$ para esos puntos) para obtener el valor de ${\displaystyle f_{xy}(x,y)}$ para el punto objetivo. O se puede hacer en la dirección opuesta, calculando primero ${\displaystyle f_{y}}$ y luego ${\displaystyle f_{x}}$ a partir de ellos, dado que ambos métodos dan resultados equivalentes.

En los bordes del conjunto de datos, cuando faltan algunos de los puntos circundantes, los puntos faltantes se pueden aproximar mediante varios métodos. Un método simple y común es asumir que la pendiente desde el punto existente hasta el punto objetivo continúa sin cambios adicionales y utilizar esto para calcular un valor hipotético para el punto faltante.

La interpolación mediante spline bicúbico requiere la solución del sistema lineal descrito anteriormente para cada celda de la cuadrícula. Se puede obtener un interpolador con propiedades similares aplicando una convolución con el siguiente núcleo en ambas dimensiones:

${\displaystyle W(x)={\begin{cases}(a+2)|x|^{3}-(a+3)|x|^{2}+1&{\text{ para }}|x|\leq 1,\\a|x|^{3}-5a|x|^{2}+8a|x|-4a&{\text{ para }}1<|x|<2,\\0&{\text{ en caso contrario}},\end{cases}}}$

donde ${\displaystyle a}$ generalmente se establece en −0,5 o −0,75. Téngase en cuenta que ${\displaystyle W(0)=1}$ y ${\displaystyle W(n)=0}$ para todos los números enteros distintos de cero ${\displaystyle n}$.

Este enfoque fue propuesto por Keys, quien demostró que ${\displaystyle a=-0.5}$ produce convergencia de tercer orden con respecto al intervalo de muestreo de la función original.^[1]​

Si se usa la notación matricial para el caso común ${\displaystyle a=-0.5}$, se puede expresar la ecuación de una manera más compacta:

${\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}}}$

para ${\displaystyle t}$ entre 0 y 1 para una dimensión. Téngase en cuenta que para la interpolación de convolución cúbica unidimensional se requieren 4 puntos de muestra. Para cada consulta se ubican dos muestras a su izquierda y dos muestras a su derecha. Estos puntos están indexados de −1 a 2 en este texto. La distancia desde el punto indexado con 0 hasta el punto de consulta se indica aquí mediante ${\displaystyle t}$.

Para dos dimensiones, se aplica primero una vez en ${\displaystyle x}$ y nuevamente en ${\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}).}$

El algoritmo bicúbico se utiliza con frecuencia para escalar imágenes y vídeos para su visualización (consúltese remuestreo de mapas de bits). Conserva los detalles finos mejor que el algoritmo de interpolación bilineal común.

Sin embargo, debido a los lóbulos negativos del núcleo, causa sobrepaso (halos). Esta circunstancia puede causar recorte y es un artefacto (véase también artefactos de anillo), pero aumenta la acutancia (nitidez aparente) y puede ser deseable.

• Portal:Matemática. Contenido relacionado con Matemática.

• Aplicación de interpolación a muestras de elevación
• Teoría de la interpolación
• Explicación e implementación Java/C++ de la interpolación (bi)cúbica
• Función de hoja de cálculo de Excel para interpolación bicúbica de Lagrange