Reescalado de imágenes con el paquete 'terra' de R

Siempre he echado en falta que R no disponga de funciones estándar para reescalar imágenes (matrices/arrays). Existen varios paquetes con esta funcionalidad, y uno de ellos que empecé a usar hoy es bueno (calidad de implementación de los algoritmos de interpolación habituales con profusión de opciones de promediado), bonito (intuitivo y fácil de usar) y barato (aceptablemente rápido).

Se trata del paquete terra, reemplazo del paquete raster, mucho menos optimizado para esta tarea. Pese a ser un paquete orientado a procesar datos GIS, su implementación del reescalado (función resample()) de objetos raster es todo lo que necesitamos para cambiar el tamaño de cualquier imagen.

~~~

Para comprobar el rendimiento del reescalado de esta función usaremos imágenes de prueba incluyendo tres tipologías, de más a menos sintéticas:

• Un test de tortura frente al aliasing: círculos concéntricos de frecuencia espacial creciente con el radio.
• Un par de imágenes digitales: mapa representado por sus contornos (de 1 píxel de ancho y sin antialiasing aplicado) y una versión "sólida" del mismo.
• Una fotografía: imagen orgánica con profusión de detalle continuo.

La primera imagen requiere algo de explicación y para visualizarla correctamente debe abrirse a tamaño real en un PC, ya que si la vemos en el móvil probablemente se reescalará de forma rápida en pantalla pudiendo aparecer moiré. Es un patrón diseñado para evaluar cómo de bueno es el filtrado paso bajo que realiza cada algoritmo, requisito imprescindible para que patrones de alta frecuencia espacial como éste no generen artefactos por aliasing tras una reducción de tamaño de la imagen (hacer clic para verla a máxima resolución):

Un correcto reescalado a la baja de esta imagen debería consistir en círculos concéntricos partiendo desde el centro, que sean sustituidos por un tono gris plano, sin artefactos ni patrones visibles, en cuanto la proximidad entre círculos resulte excesiva para poder ser codificada en la resolución de salida escogida. Cualquier artefacto que aparezca más allá de ese punto (límite de Nyquist), será una interferencia indeseable debida a un deficiente (o inexistente) filtrado paso bajo.

~~~

Las imágenes de prueba las vamos a reducir primero de sus originales 1.000x1.000 píxeles a 200x200 píxeles cada una, obteniendo el siguiente resultado ampliadas al 200% para apreciar mejor los detalles (imprescindible hacer clic para verlas en su resolución original):

Igualmente hacemos un escalado al alza, en esta ocasión reescalamos el composite de 2.000x2.000 píxeles con las cuatro imágenes a 7.547x7.547 píxeles. Seleccionamos estas zonas de interés (imprescindible hacer clic para verlas en su resolución original):

En amarillo las características más reseñables de cada algoritmo en cada tipo de imagen y reescalado. Los artefactos en los bordes de 'cubic' y 'lanczos' son difíciles de ver a simple vista, así que para resaltarlos ponemos aquí en gris medio todo lo que sea blanco o negro puros, mostrando con claridad los nuevos tonos grises que crea cada algoritmo (por eso 'near' aparece en blanco). Tanto 'cubic' como 'lanczos' generan residuos indeseables en los bordes.

Mi elección de algoritmo a priori sería:

• Para gráficos de ordenador (imágenes sintéticas): en las reducciones usaría 'bilinear' por su buen compromiso entre definición y ausencia de artefactos, pero para realizar escalados al alza me decantaría por 'cubic' por su mayor contraste, pese a que puede mostrar artefactos residuales (usando 'near' cuando necesitemos explícitamente un reescalado sin realizar interpolación, por ejemplo para no crear tonos nuevos).
• Para fotografía (lo que incluye mapas raster aéreos y satelitales) usaría siempre 'lanczos' por su mayor nivel de detalle y buena protección frente al aliasing en las reducciones. Sus artefactos no van a ser visibles en imágenes fotográficas y además tras 'cubicspline' es el algoritmo que mejor disimula el pixelado en los aumentos de tamaño.

~~~

He dejado para el final la comparativa en cuanto a velocidad de ejecución, ya que pese a existir diferencias en los tiempos todos los algoritmos funcionan aceptablemente rápido, así que difícilmente este parámetro decidirá la elección del método de interpolación.

Los tiempos obtenidos para reducciones de tamaño (desde 2.000x2.000 píxeles a 327x327 píxeles) y para escalados al alza (desde 2.000x2.000 píxeles a 7.547x7.547 píxeles) son buen proxy para intuir la complejidad de cada algoritmo:

La mejora respecto al paquete raster es de más de un orden de magnitud (unas 15 veces más rápido). Me sigue chocando no obstante que el reescalado por nearest neighbour ('near'), tenga tiempos cercanos a los verdaderos algoritmos de interpolación, cuando debería ser prácticamente instantáneo ya que solo consiste en tomar muestras de la imagen de partida sin hacer procesado.

~~~

He definido una función arrayresample() que tomando como parámetros cualquier array a reescalar, las dimensiones de salida y el algoritmo de interpolación, devuelve el array reescalado. Sirve tanto para imágenes en escala de grises (matrices) como en color (array con 3 canales de color). Repositorio con el código R: GitHub.