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.

Computación

Interpolación bicúbica en el cuadrado $[0,4]\times [0,4]$ que consta de 25 cuadrados unitarios unidos. Interpolación bicúbica según la implementación del programa Matplotlib. El color indica el valor de la función. Los puntos negros son las ubicaciones de los datos prescritos que se interpolan. Obsérvese cómo las muestras de color no son radialmente simétricas

Interpolación bilineal en el mismo conjunto de datos que el anterior. Las derivadas de la superficie no son continuas sobre los límites del cuadrado

Interpolación por el vecino más cercano en el mismo conjunto de datos que en el ejemplo anterior

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

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 $a_{ij}$ . Al hacer coincidir $p(x,y)$ con los valores de la función se obtienen cuatro ecuaciones:

$f(0,0)=p(0,0)=a_{00},$
$f(1,0)=p(1,0)=a_{00}+a_{10}+a_{20}+a_{30},$
$f(0,1)=p(0,1)=a_{00}+a_{01}+a_{02}+a_{03},$
$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 $x$ e $y$ :

$f_{x}(0,0)=p_{x}(0,0)=a_{10},$
$f_{x}(1,0)=p_{x}(1,0)=a_{10}+2a_{20}+3a_{30},$
$f_{x}(0,1)=p_{x}(0,1)=a_{10}+a_{11}+a_{12}+a_{13},$
$f_{x}(1,1)=p_{x}(1,1)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}a_{ij}i,$
$f_{y}(0,0)=p_{y}(0,0)=a_{01},$
$f_{y}(1,0)=p_{y}(1,0)=a_{01}+a_{11}+a_{21}+a_{31},$
$f_{y}(0,1)=p_{y}(0,1)=a_{01}+2a_{02}+3a_{03},$
$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 $xy$ :

$f_{xy}(0,0)=p_{xy}(0,0)=a_{11},$
$f_{xy}(1,0)=p_{xy}(1,0)=a_{11}+2a_{21}+3a_{31},$
$f_{xy}(0,1)=p_{xy}(0,1)=a_{11}+2a_{12}+3a_{13},$
$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:

p_{x}(x,y)=\textstyle \sum \limits _{i=1}^{3}\sum \limits _{j=0}^{3}a_{ij}ix^{i-1}y^{j},

p_{y}(x,y)=\textstyle \sum \limits _{i=0}^{3}\sum \limits _{j=1}^{3}a_{ij}x^{i}jy^{j-1},

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 $p(x,y)$ sobre el cuadrado unidad $[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 $a_{ij}$ en un vector

\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

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 $A\alpha =x$ .

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

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 $\alpha$ de forma rápida y sencilla.

Puede haber otra forma matricial concisa para 16 coeficientes:

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

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

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

Ampliación a retículos rectilíneos

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 $p_{x},p_{y},$ y $p_{xy}$ pasan a ser

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}},

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}},

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 $\Delta x$ es el espaciado $x$ de la celda que contiene el punto $(x,y)$ y similar para $\Delta y$ . En este caso, el enfoque más práctico para calcular los coeficientes $\alpha$ es dejar

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 $\alpha =A^{-1}x$ con $A$ como antes. A continuación, las variables de interpolación normalizadas se calculan como

{\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 $x_{0},x_{1},y_{0},$ y $y_{1}$ son las coordenadas $x$ y $y$ de los puntos de la cuadrícula que rodean el punto $(x,y)$ . Entonces, la superficie de interpolación se convierte en

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

Obtención de las derivadas de valores de funciones

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, $f_{x}$ o $f_{y}$ , usando ese método, se calcula la pendiente entre los dos puntos "circundantes" en el eje apropiado. Por ejemplo, para calcular $f_{x}$ para uno de los puntos, se debe buscar $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 $f_{y}$ .

Para encontrar la derivada cruzada $f_{xy}$ , se debe tomar la derivada en ambos ejes, uno cada vez. Por ejemplo, primero se puede usar el procedimiento $f_{x}$ para encontrar las derivadas $x$ de los puntos por encima y por debajo del punto objetivo, luego usar el procedimiento $f_{y}$ en esos valores (en lugar de, como es habitual, los valores de $f$ para esos puntos) para obtener el valor de $f_{xy}(x,y)$ para el punto objetivo. O se puede hacer en la dirección opuesta, calculando primero $f_{y}$ y luego $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.

Algoritmo de convolución bicúbica

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:

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 $a$ generalmente se establece en −0,5 o −0,75. Téngase en cuenta que $W(0)=1$ y $W(n)=0$ para todos los números enteros distintos de cero $n$ .

Este enfoque fue propuesto por Keys, quien demostró que $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 $a=-0.5$ , se puede expresar la ecuación de una manera más compacta:

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

Para dos dimensiones, se aplica primero una vez en $x$ y nuevamente en $y$ :

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

p(x,y)=p(t_{y},b_{-1},b_{0},b_{1},b_{2}).

Uso en gráficos por computadora

La mitad inferior de esta figura es una ampliación de la mitad superior, que muestra cómo se crea la aparente nitidez de la línea de la izquierda. La interpolación bicúbica provoca un sobrepaso del color, lo que aumenta la acutancia

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.

Véase también

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

Suavizado espacial
Superficie de Bézier
Interpolación bilineal
Interpolación cúbica, el análogo unidimensional del spline bicúbico
Remuestreo de Lanczos
Interpolación de vecinos naturales
Filtro Sinc
Interpolación spline
Interpolación tricúbica
Interpolación de convolución cúbica direccional

Referencias

↑ R. Keys (1981). «Cubic convolution interpolation for digital image processing». IEEE Transactions on Acoustics, Speech, and Signal Processing 29 (6): 1153-1160. Bibcode:1981ITASS..29.1153K. doi:10.1109/TASSP.1981.1163711. «citeseerx: 10.1.1.320.776».

Enlaces externos

[Keys-1] R. Keys (1981). «Cubic convolution interpolation for digital image processing». IEEE Transactions on Acoustics, Speech, and Signal Processing 29 (6): 1153-1160. Bibcode:1981ITASS..29.1153K. doi:10.1109/TASSP.1981.1163711. «citeseerx: 10.1.1.320.776».

[1]