Implementación paralela del algoritmo HSI-MSER para el registro de imágenes hiperespectrales

Implementaci´on paralela del algoritmo HSI-MSER para el registro de im´agenes hiperespectrales Daniel del Castillo1,2,´ Alvaro Ord´o˜nez3, Dora B. Heras1,2y Francisco Arg¨uello2 Resumen— El registro de im´agenes es una tarea destinada a alinear im´agenes como un paso anterior a otros tipos de procesado como son la detecci´on de cambios o el seguimiento de objetos. En el caso de im´agenes de teledetecci´on hiperespectrales tenemos una gran cantidad de informaci´on tanto espacial como espectral para cada imagen. Esto ha propiciado la aparici´on de m´etodos de registro espec´ıficos para este tipo de im´agenes que no solo usan las bandas RGB, como sucede en los algoritmos m´as cl´asicos, sino que adem´as aprovechan la informaci´on espectral contenida en todas las bandas, que pueden llegar a ser del orden de un centenar. No obstante, para algunos de estos algoritmos, especialmente los basados en detecci´on de caracter´ısticas como l´ıneas, puntos o regiones, el tiempo de ejecuci´on sigue siendo una barrera para su uso en aplicaciones que requieren un procesado especialmente r´apido. Uno de estos algoritmos es Hyperspectral Image-Maximally Stable Extremal Regions (HSI-MSER), que trata de encontrar regiones comunes en las im´agenes explotando la informaci´on disponible en las bandas espectrales para conseguir un mejor registro. En este art´ıculo se presenta una versi´on paralela del algoritmo HSI-MSER sobre una arquitectura heterog´enea de bajo coste. Ha sido dise˜nada para explotar eficientemente las arquitecturas de la CPU y GPU usando OpenMP y CUDA, respectivamente. Como resultado, para im´agenes hiperespectrales de teledetecci´on disponibles en la bibliograf´ıa estado del arte se consiguen aceleraciones de hasta 7×respecto a la versi´on secuencial de HSI-MSER. Palabras clave— Registro de im´agenes, hiperespectral, teledetecci´on, CUDA, GPU, OpenMP I. Introducci´ on EL registro de im´agenes es una operaci´on crucial a la hora de trabajar con im´agenes de teledetecci´on de superficie terrestre capturadas por drones o sat´elites. Se puede definir como un proceso que toma dos im´agenes que capturan informaci´on de la misma ´area geogr´afica y modifica la escala, la rotaci´on y el desplazamiento de una de ellas para que se alinee con la otra. Esta t´ecnica puede ser necesaria para diferentes prop´ositos. El m´as habitual es realizar un an´alisis de cambios sobre im´agenes de la misma ´area capturadas en distintos momentos, a diferentes alturas o con equipos distintos. Es importante destacar que existe una gran variedad de algoritmos de registro. La mayor´ıa de ellos 1Centro Singular de Investigaci´on en Tecnolox´ıas Intelixentes (CiTIUS), Universidade de Santiago de Compostela, email: {d.delcastillo,dora.blanco}@usc.gal. 2Departamento de Electr´onica e Computaci´on, Universidade de Santiago de Compostela, e-mail: [email protected]. 3Universidade da Coru˜na, Grupo Integrado de Ingenier´ıa, CITIC, Elvi˜na, 15071 A Coru˜na, e-mail: [email protected] se pueden clasificar en dos tipos: basados en ´area o basados en caracter´ısticas. Los m´etodos basados en ´area, como los basados en la Fourier-Mellin transform (FMT) [1], trabajan directamente con las intensidades buscando correlaciones para realizar el alineamiento. Por el contrario, los m´etodos basados en caracter´ısticas, como Speeded Up Robust Features (SURF) [2], buscan ciertos elementos dentro de cada imagen, como l´ıneas o puntos. Esto hace que tengan un coste computacional mayor que los basados en intensidad, pero, sin embargo, consiguen una mayor tolerancia a cambios en la imagen [3]. Esto permite registrar im´agenes con mayores diferencias de escala entre ellas, lo que resulta de gran utilidad para ciertas aplicaciones cuando se ha de trabajar con im´agenes obtenidas, por ejemplo, a alturas muy diferentes. En la actualidad, cada vez son m´as los sistemas que incorporan sensores hiperespectrales que permiten obtener mapas detallados del terreno con una alta resoluci´on espectral, lo que hace necesario disponer de t´ecnicas de registro que hagan un uso eficiente y eficaz de toda la informaci´on que contienen. HSI-MSER [4] es un algoritmo para el registro de dos im´agenes hiperespectrales basado en la detecci´on de caracter´ısticas, en este caso, regiones de inter´es. El m´etodo explota eficientemente la informaci´on espectral disponible en distintas bandas espectrales e incorpora este tipo de informaci´on en el descriptor que se construye para cada regi´on detectada. Esto permite registrar im´agenes con diferencias de escala de hasta 15×. A pesar de ser m´as tolerantes a cambios, los m´etodos basados en caracter´ısticas presentan como desventaja un mayor tiempo de ejecuci´on, ya que se basan en una extracci´on profunda de caracter´ısticas como, por ejemplo, regiones con determinadas propiedades sobre ambas im´agenes. Las caracter´ısticas son a continuaci´on descritas mediante vectores de gran longitud para realizar despu´es un proceso de emparejamiento entre las dos im´agenes que requiere una b´usqueda exhaustiva. Esta desventaja se hace notar especialmente a la hora de trabajar con im´agenes hiperespectrales, ya que una primera aproximaci´on exigir´ıa buscar caracter´ısticas sobre todas las bandas espectrales de ambas im´agenes. Sin embargo, las empresas y usuarios que hacen uso de este tipo de im´agenes no siempre tienen acceso a supercomputadores o infraestructuras computacionales costosas. Adem´as, resulta mucho m´as conveniente en muchos Jornadas SARTECO 2023 331 casos tratar este tipo de im´agenes directamente en el punto de recogida de los datos para poder ofrecer a los expertos informaci´on que facilite una toma de decisiones r´apida. Es por ello que se necesitan algoritmos especialmente dise˜nados para explotar el paralelismo del hardware com´unmente disponible sobre el terreno. En la actualidad este podr´ıa, por ejemplo, consistir en un sistema multin´ucleo provisto de una o varias GPUs. Existen en la literatura diferentes algoritmos dise˜nados para su ejecuci´on en este tipo de arquitecturas paralelas, pero ninguno de ellos est´a basado en registro mediante detecci´on de regiones de inter´es [5], [6]. En este art´ıculo se presenta una primera implementaci´on paralela del algoritmo HSI-MSER [4], que realiza registro de im´agenes hiperespectrales basado en la detecci´on de regiones delimitadas mediante elipses, para su ejecuci´on en un sistema heterog´eneo basado en un procesador multicore y GPU. El algoritmo ha sido desarrollado en OpenMP y CUDA para obtener una explotaci´on eficiente de la arquitectura. El resto del art´ıculo tiene la siguiente estructura. La secci´on II explica las bases de este trabajo y el algoritmo original HSI-MSER. La secci´on III describe la implementaci´on h´ıbrida en CPU y GPU desarrollada. En la secci´on IV se pueden encontrar una evaluaci´on de los resultados conseguidos, analizando no solamente la eficiencia en cuanto a ahorro de tiempo de computaci´on, sino tambi´en la calidad en cuanto a capacidad de registro de los algoritmos paralelos desarrollados. Por ´ultimo, en la secci´on V se indican los retos alcanzados, as´ı como v´ıas de trabajo futuro. II. Algoritmo original En esta secci´on se describe el algoritmo HSI-MSER para el registro de dos im´agenes hiperespectrales de teledetecci´on, para luego explicar el dise˜no del algoritmo paralelo mediante OpenMP y CUDA. Maximally stable extremal regions (MSER) [7] es un algoritmo de detecci´on de caracter´ısticas que extrae regiones que persisten al analizar la imagen con diferentes umbrales de gris. Estas regiones son invariantes a las transformaciones de la imagen, lo que las hace ideales para el registro. Scale Invariant Feature Transform (SIFT) [8] es un conocido detector y descriptor de puntos clave basado en la construcci´on de un espacio de escala gaussiano. HSI-MSER adapta MSER para su uso en im´agenes hiperespectrales, utilizando una adaptaci´on de SIFT como descriptor espacial de regiones e introduciendo un m´etodo para seleccionar qu´e bandas espectrales utilizar en funci´on de su entrop´ıa, denominado Entropy-Based Band Selection (EBS). HSI-MSER adem´as introduce un descriptor espectral a cada regi´on detectada y propone un m´etodo exhaustivo para el c´alculo de la transformaci´on geom´etrica que al ser aplicada sobre una imagen la alinea con la otra. Todo esto ha permitido a este algoritmo aprovechar la informaci´on presente en las distintas bandas de una imagen hiperespectral para mejorar la precisi´on de registro. El algoritmo HSI-MSER est´a compuesto por 6 fases, como se puede ver en la Figura 1, que son explicadas a continuaci´on. Selecci´on de bandas. Las im´agenes hiperespectrales est´an formadas por cientos de bandas contiguas que, dada su proximidad en cuanto a la longitud de onda asociada a ellas, contienen informaci´on redundante en algunos casos. Habitualmente se aplica un m´etodo de extracci´on de las bandas estad´ısticamente m´as relevantes que permita realizar el registro sobre un n´umero de ellas que no supere la decena. Tal como se puede observar en la Figura 1, la primera fase del algoritmo consiste en aplicar el m´etodo de selecci´on de bandas llamado Entropy-Based Band Selection (EBS) [4]. Este m´etodo tiene la peculiaridad de tener en cuenta las dos im´agenes, en lugar de escoger las bandas de cada imagen por separado. EBS selecciona las Nbandas con mayor entrop´ıa que est´an separadas por al menos Dbandas, cuando estas est´an ordenadas por longitud de onda. De esta forma consigue seleccionar bandas con una alta entrop´ıa, pero que tambi´en difieren sustancialmente en longitud de onda. Detecci´on de regiones. En esta segunda fase, se analizan cada una de las bandas escogidas en ambas im´agenes en busca de regiones de inter´es. Para ello se usa una adaptaci´on del algoritmo MSER [7]. Este algoritmo procesa todos los p´ıxeles de la imagen en orden de intensidad y busca regiones que mantengan su forma al ir aumentando el umbral de nivel de gris. De esta forma se obtienen regiones que existen en distintas bandas espectrales. Despu´es, estas regiones son definidas como elipses. Descripci´on de las regiones. En esta fase se describen las regiones que se detectaron durante la fase anterior. A cada regi´on se le asigna un descriptor, es decir, una secuencia de n´umeros que resumir´a distintas caracter´ısticas de la misma. Las mismas regiones, pero extra´ıdas en diferentes im´agenes y bajo cambios de iluminaci´on o transformaciones geom´etricas distintas, deben dar lugar a un descriptor similar. Esto nos permitir´a emparejar las regiones de ambas im´agenes en la siguiente fase. HSI-MSER propone un descriptor compuesto por informaci´on tanto espacial como espectral. El algoritmo SIFT [8] es usado para construir la parte espacial del descriptor. Como SIFT fue dise˜nado para describir puntos y no regiones, en HSIMSER se emplea una adaptaci´on que tiene en cuenta el ´area de cada regi´on. Para computar el descriptor primero se calculan los ´angulos dominantes de cada regi´on usando el gradiente de su superficie y una ventana circular con ponderaci´on gaussiana. Una vez se ha determinado el ´angulo dominante se puede construir el descriptor, compuesto por 128 par´ametros. Estos valores se normalizan para reducir la influencia de cambios en la iluminaci´on o el brillo. Adem´as de este descriptor espacial, tambi´en se tiene en cuenta la informaci´on espectral de cada regi´on. Para ello se usa la firma espectral del centro de la elipse de todas las bandas, no solo de las escogidas durante la selecci´on de bandas. 332 Jornadas SARTECO 2023 Fig. 1: Plataforma de computaci´on paralela utilizada en cada una de las fases del algoritmo HSI-MSER. Adaptado de [4]. Emparejamiento de regiones. Una vez que las regiones han sido descritas, tratamos de encontrar parejas entre las regiones de ambas im´agenes. Estos emparejamientos se hacen de forma independiente para cada banda de la imagen, pero se tiene en cuenta la informaci´on espectral del resto de bandas, como se explic´o en la subsecci´on anterior. Para calcular la similitud entre dos regiones que determina si dos regiones est´an emparejadas se usa la distancia eucl´ıdea para la parte espacial y la similitud coseno para la parte espectral. ´ Unicamente se considera un emparejamiento si estas distancias est´an dentro de unos par´ametros establecidos [4]. Combinaci´on de las bandas. En este paso se combinan todas las parejas de regiones provenientes de todas las bandas. Esto permite usar en el registro regiones que solo aparecen en algunas frecuencias del espectro. Registro. Para registrar las im´agenes, HSI-MSER realiza una b´usqueda exhaustiva basada en histogramas [4]. Para cada pareja se computan el ´angulo, la escala y el desplazamiento correspondientes, tal y como se puede ver en la Figura 1. Los distintos ´angulos se guardan en un histograma que representa los 360°, divididos en 72 contenedores con un tama˜no base de 5°y un solapamiento de 2.5°a cada lado, llevando el total a 10°de tama˜no. Una vez que este histograma se ha construido, las escalas del contenedor con m´as elementos se ordenan para obtener la pareja o parejas que representan la mediana. Estas parejas de regiones son las que se emplean para calcular la transformaci´on geom´etrica que nos permite alinear las im´agenes hiperespectrales. III. Algoritmo paralelo propuesto para registro mediante HSI-MSER En esta secci´on, se presenta la implementaci´on paralela de HSI-MSER en CUDA y OpenMP para explotar arquitecturas heterog´eneas. Se ha escogido para cada fase del algoritmo el paradigma de computaci´on m´as adecuado con el fin de maximizar rendimiento. En el Algoritmo 1 se presenta el pseudoc´odigo de la implementaci´on paralela de HSI-MSER. Para cada fase, se indica con <y>que partes han sido paralelizadas en CUDA y con # las partes que lo han sido con OpenMP. A continuaci´on se explican los detalles asociados a cada fase. Selecci´on de bandas. Esta fase se paraleliz´o en la GPU usando el modelo de programaci´on CUDA. La naturaleza completamente paralela del c´alculo de entrop´ıa hace que este proceso sea realizable en la arquitectura de una GPU de manera eficiente [9]. Para el c´omputo de la entrop´ıa de cada banda (l´ınea 1, Alg. 1) se usan varias funciones de la biblioteca CUB [10]. En concreto, DeviceReduce::Min yDeviceReduce::Max se utilizan para calcular de forma eficiente los valores m´aximo y m´ınimo de cada banda. Una vez se tienen estos valores, se usa DeviceHistogram::HistogramEven para calcular el histograma. Por ´ultimo, se realiza una reducci´on de los valores del histograma, haciendo uso de las operaciones a nivel de warp para mejorar su rendimiento. Cuando ya se tienen los valores para la entrop´ıa de cada banda, se procede a seleccionar las bandas en CPU (l´ınea 2, Alg. 1). Para ello se ordenan las bandas y se procesan siguiendo lo explicado en la Subsecci´on II. Extracci´on de regiones. MSER se basa en el algoritmo Union find. Este algoritmo tiene una complejidad de O(m·α(m, n)), donde nes el n´umero de elementos, mes el n´umero de operaciones de b´usqueda que se necesitar´an y αes la funci´on inversa de Ackerman [11]. mes mayor o igual que n, pero tiene una relaci´on de car´acter lineal (m=kn). Esta complejidad se considera como cuasi lineal debido a la lentitud con la que crece la funci´on inversa de Ackerman. En la literatura ya se han publicado t´ecnicas que se pueden aplicar para paralelizar este algoritmo [12], incluso en la GPU [13], [14]. Sin embargo, MSER, a pesar de estar basado en Union find, presenta varias diferencias respecto a este. Una de ellas es que los p´ıxeles deben ser procesados por orden de intensidad creciente, lo que dificulta la paralelizaci´on de la construcci´on del ´arbol, ya que cada p´ıxel es colocado en el ´arbol respecto a las posiciones de los p´ıxeles anteriores. Si se trata de dividir la imagen para realizar la construcci´on del ´arbol en paralelo, ya sea por posici´on f´ısica de cada p´ıxel o por intensidad, se consiguen ´arboles inconexos. En [15] se propone una Jornadas SARTECO 2023 333 Algorithm 1 Pseudoc´odigo de la implementaci´on h´ıbrida en OpenMP y CUDA Entrada: Imagen de referencia I1e imagen objetivo I2. Salida: Escala ρ, ´angulo θy desplazamiento (x,y). Fase I. Selecci´on de bandas. 1: <C´alculo de la entrop´ıa para cada banda > 2: Selecci´on de las bandas con base en la entrop´ıa 3: para cada banda bseleccionada: Fase II. Detecci´on de regiones. 4: Realizar particiones para cada imagen seg´un el n´umero de hilos disponibles # Paralelizado con OpenMP 5: Detectar regiones por separado en cada partici´on Fase III. Descripci´on de regiones. 6: <Calcular la orientaci´on dominante de cada regi´on > 7: <Computar una ventana gaussiana a partir de la banda > 8: <Computar el descriptor de cada regi´on > Fase IV. Emparejamiento de regiones. 9: <Calcular las distancias entre los descriptores espaciales de las regiones de ambas im´agenes > 10: para cada regi´on rde Ib 1: 11: Si La menor distancia desde ra una regi´on de Ib 2cumple ciertos par´ametros 12: Si La similitud espectral cumple ciertos par´ametros 13: Guardar la pareja Fase V. Combinaci´on de bandas. 14: Combinar las parejas encontradas en distintas bandas en un ´unico vector Fase VI. Registro. 15: Computar una transformaci´on geom´etrica para cada una de las parejas 16: Generar un histograma de ´angulos 17: Escoger el rango de ´angulos m´as repetido y usar la pareja mediana para el registro final Fig. 2: Esquema de paralelizaci´on de la fase de extracci´on de regiones. paralelizaci´on de MSER que fusiona estos ´arboles resultantes, pero que no se almacena el ´area de cada regi´on, dato necesario para convertir estas regiones en elipses. La Figura 2 muestra en detalle c´omo ha sido finalmente paralelizada esta etapa. Primero, se divide la imagen horizontalmente en un n´umero de particiones igual al n´umero de hilos disponibles (l´ınea 4, Alg. 1). Luego, como se puede ver en la Figura 2, se buscan las regiones en cada partici´on concurrentemente (l´ınea 5, Alg. 1) haciendo uso de OpenMP [16]. Como la detecci´on de regiones no modifica la imagen, la memoria puede ser compartida y no hay necesidad de comunicaci´on expl´ıcita entre hilos. Cada hilo 334 Jornadas SARTECO 2023 de c´omputo guarda localmente las regiones que encuentra, hasta que, finalmente, se unen todas. Para evitar que no se detecten incorrectamente regiones que est´en situadas justo en los l´ımites de una partici´on, se ha a˜nadido un solapamiento entre dichas particiones. Este solapamiento se ha ilustrado en la Figura 2 con el color rojo. Un 20 % de solapamiento significa que cada hilo procesa un 120 % respecto a la versi´on sin solape, compartiendo un 10 % con cada uno de los hilos vecinos. Debido a la extracci´on de regiones por particiones, ha sido necesario ajustar el par´ametro de min area de HSI-MSER para que tenga en cuenta el tama˜no de la imagen completa y no solo el de la partici´on. Este par´ametro determina cu´al es el tama˜no m´ınimo de una regi´on para ser considerada. Su valor representa un porcentaje respecto al tama˜no total de la imagen, es decir, si el valor es un 10 %, cada regi´on tiene que tener al menos el 10 % de la imagen para ser considerada. La idea detr´as de este c´alculo es que, im´agenes de menor tama˜no es probable que tengan regiones m´as peque˜nas. El par´ametro min area ha sido ajustado para que detecte regiones del mismo tama˜no que si se procesara la banda completa. Descripci´on de las regiones. El proceso de descripci´on de cada regi´on se ha paralelizado en la GPU, tal y como aparece en la Figura 2. Para ello se han implementado tres kernels CUDA. En concreto, el c´ompu- to del histograma de los ´angulos (l´ınea 6, Alg. 1), la generaci´on de la ventana de ponderaci´on gausiana (l´ınea 7, Alg. 1) y el c´omputo final de los descriptores (l´ınea 8, Alg. 1). Estos kernels se caracterizan por una alta carga de operaciones aritm´eticas y trigonom´etricas. Adem´as, excepto en el caso de la ventana gaussiana, es necesario el uso de instrucciones at´omicas. Estos kernels han sido analizados y optimizados siguiendo la metodolog´ıa Assess, Parallelize, Optimize, Deploy (APOD) [17] y haciendo uso de la herramienta de profiling NVIDIA Nsight Compute. Se han aplicado las optimizaciones recomendadas por NVIDIA [17] para obtener el mejor rendimiento en la GPU, por ejemplo, minimizar la transferencia de datos entre la memoria principal y GPU, escoger el n´umero de hilos por bloque y el n´umero de bloques ´optimo por cada kernel, y maximizar la ocupancia de los multiprocesadores evitando divergencias dentro de los warps de cada kernel, entre otras. Emparejamiento de regiones. Para el emparejamiento se ha usado una aproximaci´on de la distancia eucl´ıdea [9]. Este m´etodo traduce el c´alculo de las distancias entre las regiones de una misma banda (l´ınea 9, Alg. 1) en operaciones matriciales. Las matrices son ideales para su c´omputo en GPU, ya que se procesan repitiendo la misma operaci´on y conforman un conjunto de datos sin dependencias. Esto nos permite obtener una gran aceleraci´on en esta etapa. Adem´as, se ha empleado la biblioteca cuBLAS [18] siguiendo la implementaci´on descrita en [9]. Una vez se tienen las distancias entre los descriptores espaciales, se eval´ua la distancia espectral para aquellos que cumplan ciertos par´ametros establecidos por el algoritmo HSI-MSER original y se guarda la pareja si la similitud espectral tambi´en es satisfactoria (l´ıneas 10-13, Alg. 1). Combinaci´on de las bandas y registro. Ambas fases se ejecutan secuencialmente en la CPU debido a su bajo coste computacional, tal y como se mostrar´a en la siguiente secci´on (l´ıneas 14-17, Alg. 1). IV. Resultados En esta secci´on se presentan las condiciones de experimentaci´on y las im´agenes utilizadas, as´ı como los resultados obtenidos tanto en t´erminos de tiempos de ejecuci´on y rendimiento como en t´erminos de capacidad de registro. El algoritmo ha sido ejecutado sobre un ordenador con un procesador Intel Core i7 9700k de 3.6 GHz con 8 cores, 32 GB de RAM y una GPU NVIDIA GeForce RTX 3050. Por otro lado, se han usado g++ 11.3 y nvcc 12.1 como compiladores y Ubuntu 22.04.2 LTS como sistema operativo. Las versiones de las bibliotecas de NVIDIA como CUB y cuBLAS vienen determinadas por la versi´on 12.1 de CUDA. Los tiempos de ejecuci´on y aceleraciones son el resultado de promediar 10 ejecuciones independientes. Por ´ultimo, se han usado las opciones O3 2 arch=native”para generar c´odigo optimizado para la arquitectura correspondiente. El algoritmo HSI-MSER se ha ejecutado con los valores por defecto [4] a excepci´on del par´ametro min area tal y como se explic´o en la secci´on III. Se han empleado 14 im´agenes para realizar los experimentos de este trabajo. Las caracter´ısticas de cada una de ellas pueden encontrarse en el art´ıculo del m´etodo original [4] y est´an disponibles para su descarga en [19]. Las im´agenes Pavia University yPavia Centre han sido tomadas por el sensor Reflective Optics System Imaging Spectrometer (ROSIS). Las restantes han sido capturadas por el Airborne Visible /Infrared Imaging Spectrometer (AVIRIS) [20]. Para Jasper Ridge,Santa Barbara Front,Santa Barbara Line,Baraboo Hills yCrown Point se cuenta con dos im´agenes que han sido capturadas en distintos a˜nos y, por lo tanto, presentan cambios de escala, rotaci´on, desplazamiento e incluso distorsi´on. Nos referiremos a estas im´agenes como casos reales. La escala, ´angulo y desplazamientos de referencia entre estas im´agenes est´an disponibles en [4]. Para trabajar con las im´agenes sin pareja, se genera una imagen objetivo aplic´andole una escala y un ´angulo a la original, lo que nos permite estudiar el registro bajo condiciones controladas. A. Tiempos de ejecuci´on y aceleraciones La Tabla I presenta los tiempos de ejecuci´on de la implementaci´on original secuencial y de la implementaci´on paralela propuesta en este trabajo utilizando distinto n´umero de hilos para los casos reales. La ver- si´on CUDA ejecuta la selecci´on de bandas, la descripci´on de regiones y el emparejamiento en la GPU, lo que le da una ventaja considerable frente a la versi´on secuencial. Las versiones CUDA + OpenMP, adem´as Jornadas SARTECO 2023 335 Tabla I: Tiempos de ejecuci´on (segundos) de la implementaci´on secuencial y de la implementaci´on paralela propuesta con distinto n´umero de hilos en los casos reales. Imagen Secuencial CUDA CUDA + CUDA + CUDA + OpenMP 2 OpenMP 4 OpenMP 8 Santa Barbara Front 38.24 7.22 5.54 4.64 4.27 Crown Point 39.96 11.38 7.05 5.78 5.21 Jasper Ridge 50.41 12.35 7.79 5.98 5.33 Santa Barbara Line 36.69 12.00 7.65 5.91 5.14 Baraboo Hills 36.07 12.00 7.32 5.91 4.99 Media 40.27 10.45 6.26 4.86 4.22 de ejecutar las partes mencionadas en la GPU, ejecutan la detecci´on de regiones usando 2, 4 y 8 hilos con la ayuda de OpenMP. La Figura 3 muestra las aceleraciones de las implementaciones paralelas frente la implementaci´on secuencial. Los resultados muestran que gracias a la paralelizaci´on de la selecci´on de bandas, la descripci´on de regiones y el emparejamiento en GPU, el tiempo medio se reduce de un valor promedio de 40.27s a 10.45s. Esto hace que se consiga una aceleraci´on promedio de 4×. Si adem´as a˜nadimos la paralelizaci´on de la detecci´on de regiones en OpenMP, en el caso con 8 hilos, se consigue una aceleraci´on cercana a 10×. Para esta versi´on, el tiempo promedio para registrar dos im´agenes hiperespectrales es de 4.22s, lo que constituye una mejora notable respecto a los 40.27s del m´etodo secuencial. Fig. 3: Aceleraciones de la implementaci´on paralela respecto a la secuencial en distintas im´agenes. Se utilizan los tiempos de la Tabla I. B. Capacidad de registro Para comprobar c´omo ha afectado la paralelizaci´on a la capacidad de registro, se han llevado a cabo varios experimentos. Se ha seguido el mismo procedimiento propuesto en el art´ıculo original de HSIMSER [4]. Este m´etodo se basa en ir aumentando y disminuyendo la escala original de cada imagen a la vez que se aplican rotaciones con el fin de aumentar el conjunto de im´agenes de prueba. La escala se va aumentando en incrementos de 0.5×y se reduce en fracciones (1/2, 1/3, 1/4...). Para cada escala se aplica una rotaci´on de 0 a 355°, en incrementos de 0.5°. Es decir, para cada incremento de escala estamos probando 72 rotaciones. Se considera que una escala ha sido correctamente registrada cuando se registran todos los ´angulos asociados a la misma. Esto hace un total de 2592 casos probados por imagen. La Tabla II resume el rango de escalas registradas con ´exito para cada una de las im´agenes utilizando las distintas versiones del algoritmo. Se indica entre par´entesis el n´umero de escalas correctamente registradas para todos los 72 ´angulos. Como se puede observar, el n´umero de escalas correctamente registradas disminuye al aumentar el n´umero particiones. El algoritmo secuencial es capaz de registrar una media de 19 escalas en comparaci´on con las 10.89 de la ver- si´on m´as eficiente en t´erminos de c´omputo (CUDA + OpenMP 8). Al aumentar el n´umero de hilos, aumenta el n´umero de particiones horizontales en las que se dividen las bandas en la fase de extracci´on de regiones. Esta divisi´on de datos entre hilos deteriora la capacidad de detectar regiones en los bordes, reduciendo el n´umero de emparejamientos correctos y, por tanto, impidiendo realizar un registro adecuado en casos con una gran diferencia de escala. Habiendo estudiado los tiempos de ejecuci´on, aceleraciones y el n´umero de escalas correctamente registradas en comparaci´on la implementaci´on secuencial original, se selecciona la implementaci´on CUDA + OpenMP 4 como una soluci´on de compromiso entre aceleraci´on y capacidad de registro. C. Influencia de la adici´on de solapamiento de datos entre particiones sobre la calidad de registro Como se explic´o en la Secci´on III, cada banda de la imagen se divide horizontalmente en tantas particiones como n´umero de hilos disponibles. Cada hilo extrae regiones de forma paralela en una partici´on distinta. Para evitar la p´erdida de regiones entre las particiones, efecto analizado en la secci´on anterior, se a˜nade un solapamiento, es decir, una regi´on frontera a˜nadida a los datos que se asignan a cada hilo. La Tabla III muestra el n´umero de escalas correctamente registradas para distintos tama˜nos de solapamiento (0 %, 10 %, 20 % y 30 % respecto al tama˜no de la partici´on asignada a cada hilo) con la implementaci´on CUDA + OpenMP 4, que fue seleccionada como decisi´on de compromiso. Se puede observar que el n´umero de escalas correctamente registradas incrementa al aumentar el solapamiento. Sin embargo, podemos ver que un solapamiento del 30 % obtiene los mismos resultados que el de 20 %, una media de 15.44 escalas en comparaci´on con las 13.44 obte- 336 Jornadas SARTECO 2023 Tabla II: Capacidad de registro. Rango y n´umero de escalas registradas con ´exito para cada imagen para las distintas implementaciones del algoritmo. Imagen Secuencial [4] CUDA + OpenMP 2 CUDA + OpenMP 4 CUDA + OpenMP 8 Pavia University 1/7×- 12.0×(29) 1/5×- 10.5×(24) 1/5×- 9.0×(21) 1/5×- 7.5×(18) Pavia Centre 1/8×- 15.0×(36) 1/8×- 13.0×(32) 1/7×- 10.0×(25) 1/6×- 7.5×(19) Indian Pines 1/3×- 8.0×(17) 1/3×- 6.5×(14) 1/2×- 5.0×(10) 1/1×- 4.0×( 7) Salinas 1/5×- 7.0×(17) 1/5×- 6.5×(16) 1/4×- 7.0×(16) 1/3×- 5.5×(12) Jasper Ridge 1/3×- 7.0×(15) 1/2×- 3.5×( 7) 1/2×- 3.0×( 6) 1/2×- 3.0×( 6) Santa Barbara Front 1/3×- 7.0×(15) 1/2×- 7.0×(14) 1/2×- 5.5×(11) 1/2×- 5.5×(11) Santa Barbara Line 1/5×- 7.0×(17) 1/4×- 6.0×(14) 1/4×- 5.0×(12) 1/3×- 4.0×( 9) Baraboo Hills 1/2×- 4.0×( 8) 1/1×- 3.5×( 6) 1/1×- 4.0×( 7) 1/1×- 3.0×( 5) Crown Point 1/7×- 6.0×(17) 1/5×- 6.0×(15) 1/4×- 5.5×(13) 1/4×- 4.5×(11) Media 19.00 15.78 13.44 10.89 (a) Sin solapamiento. (b) 20 % de solapamiento. Fig. 4: Influencia de la adici´on de solapamiento en la detecci´on de regiones en Pavia University con 4 particiones. nidas sin solapamiento. Es por ello que se selecciona la implementaci´on CUDA + OpenMP 4 con 20 % de solapamiento. Los efectos del solapamiento se muestran en la Figura 4. En esta figura se han representado varias parejas de regiones sobre la imagen de Pavia University con una diferencia de escala de 1.75×. El ´area abarcada por el solapamiento est´a marcada en amarillo, en este caso, un 20 %. En la figura se representa c´omo ciertas regiones, pintadas en color verde, no ser´ıan detectadas sin la adici´on de solapamiento al quedar divididas por las separaciones entre particiones. D. N´umero de parejas y regiones En la secci´on anterior hemos visto como el solapamiento nos permite mejorar el n´umero de escalas correctamente registradas, haciendo que la implementaci´on paralela obtenga resultados m´as pr´oximos a la implementaci´on secuencial en cuanto a capacidad de registro de escalas. En esta secci´on se analizan los motivos de esta mejora. En la Tabla IV se comparan la implementaci´on secuencial y las implementaciones seleccionadas, la implementaci´on CUDA + OpenMP 4 con y sin solapamiento del 20 %, en base a diferentes m´etricas. Para esta comparaci´on se ha registrado cada caso real sin aplicar ning´un escalado o rotaci´on adicional. En la primera fila de cada imagen se muestra el n´umero de parejas de regiones encontradas durante la fase de emparejamiento. Se puede observar, que el n´umero de parejas disminuye entre la versi´on secuencial y la versi´on con cuatro particiones. Esto es esperable, ya que como se explic´o anteriormente, las particiones pueden dividir regiones en dos, dando lugar a su p´erdida o dos regiones de menor tama˜no, que luego pueden ser filtradas por el propio algoritmo o pueden no ser fiables. Sin embargo, se puede ver que con un 20 % de solapamiento se consigue recuperar gran parte del n´umero de parejas de la versi´on secuencial. Este mismo efecto lo vemos en la segunda y tercera fila, donde se muestran el n´umero de regiones extra´ıdas para cada imagen. El mayor n´umero de regiones obtenido en algunas im´agenes se debe a que la misma o similar regi´on es obtenida en las diferentes particiones debido al solapamiento. En la cuarta fila, se muestra el porcentaje de parejas correctas. Para calcular esta m´etrica, se mide la distancia entre las regiones que conforman cada pareja tras aplicarle a la regi´on de la imagen a registrar la transformaci´on geom´etrica de referencia [4]. De esta manera, si la distancia tras aplicar la escala, ´angulo y desplazamientos de referencia es mayor que 2 p´ıxeles, la pareja es considerada como incorrecta. Se puede ver que el porcentaje de parejas correctas Jornadas SARTECO 2023 337 Tabla III: Capacidad de registro. Rango y n´umero de escalas registradas con ´exito para cada imagen para distintos valores de solapamiento en la versi´on CUDA + OpenMP 4. Imagen 0 % 10 % 20 % 30 % Pavia University 1/5×- 9.0×(21) 1/4×- 10.0×(22) 1/4×- 11.0×(24) 1/4×- 11.0×(24) Pavia Centre 1/7×- 10.0×(25) 1/7×- 10.5×(26) 1/6×- 11.5×(27) 1/6×- 11.5×(27) Indian Pines 1/2×- 5.0×(10) 1/2×- 5.0×(10) 1/2×- 6.5×(13) 1/2×- 6.5×(13) Salinas 1/4×- 7.0×(16) 1/4×- 7.0×(16) 1/4×- 7.0×(16) 1/4×- 7.0×(16) Jasper Ridge 1/2×- 3.0×( 6) 1/2×- 5.0×(10) 1/2×- 5.0×(10) 1/2×- 5.0×(10) Santa Barbara Front 1/2×- 5.5×(11) 1/2×- 7.0×(14) 1/3×- 7.0×(15) 1/3×- 7.0×(15) Santa Barbara Line 1/4×- 5.0×(12) 1/4×- 6.0×(14) 1/4×- 6.5×(15) 1/4×- 6.5×(15) Baraboo Hills 1/1×- 4.0×( 7) 1/1×- 4.0×( 7) 1/1×- 4.0×( 7) 1/1×- 4.0×( 7) Crown Point 1/4×- 5.5×(13) 1/3×- 6.0×(13) 1/3×- 5.5×(12) 1/3×- 5.5×(12) Media 13.44 14.67 15.44 15.44 est´a, por lo general, en torno al 75 %. Este porcentaje mejora ligeramente al a˜nadir solapamiento, indicando que se recuperan m´as regiones correctas que incorrectas con su adici´on. En el caso de Crown Point, el porcentaje de parejas correctas es considerablemente menor debido a la distorsi´on no lineal que presenta. Para poder evitar esto, se necesitan m´etodos capaces de calcular transformaciones geom´etricas con mayores grados de libertad que escala, ´angulo y desplazamiento. En la quinta fila se indica el tiempo de ejecuci´on total (media de 10 ejecuciones del registro). Por ´ultimo, se ha a˜nadido el n´umero de escalas que se registran para cada uno de estos casos reales si se var´ıan escalas y ´angulos de manera exhaustiva. Al comparar las distintas implementaciones, se puede observar que existe una relaci´on clara entre el n´umero de parejas detectadas y el n´umero de escalas registradas. Cuando la diferencia de escala es demasiado grande, el registro falla por no contar con parejas para calcular la transformaci´on o porque las pocas parejas obtenidas son incorrectas. Gracias al solapamiento del 20 %, se recuperan regiones y parejas que se encontraban en los bordes de las particiones, consiguiendo aumentar el n´umero de escalas registradas. No obstante, la versi´on con solapamiento requiere de un mayor tiempo de ejecuci´on, ya que al a˜nadirle un solapamiento a cada partici´on, estamos aumentando el espacio donde debemos buscar regiones. Por consiguiente, tal y como se ha explicado, se obtiene un mayor n´umero de regiones, que deriva en un aumento del tiempo de c´omputo de las siguientes etapas (descripci´on, emparejamiento, combinaci´on y registro). En el caso m´as desfavorable, para las im´agenes Jasper Ridge, se necesitan 5.72s frente a los 5.22s de la versi´on sin solapamiento, pasando de registrar 10 a 6 escalas, respectivamente. Si lo comparamos con la versi´on secuencial, esta necesita 43.95s. E. Aceleraci´on por etapas Las Tablas V y VI presentan los tiempos de ejecuci´on y aceleraciones para las implementaciones secuencial y la seleccionada (CUDA + OpenMP 4 con 20 % de solapamiento) para el caso real de mayor tama˜no (Baraboo Hills) y para el de menor tama˜no (Santa Barbara Front). Como se puede observar, las fases m´as costosas son la detecci´on, descripci´on y emparejamiento de regiones. En ellas es donde se ha conseguido mejores aceleraciones, destacando especialmente la fase de emparejamiento que consigue una aceleraci´on pr´oxima a 21×. La aproximaci´on de la distancia eucl´ıdea en c´alculos matriciales para el emparejamiento de las regiones, nos permite explotar adecuadamente el modelo de paralelismo y arquitectura de la GPU. Las otras dos fases que se ejecutan en GPU, selecci´on de bandas y descripci´on de regiones, han obtenido aceleraciones de hasta 2.68×y 16.19×, respectivamente. La fase de detecci´on de regiones est´a paralelizado usando OpenMP y consigue una aceleraci´on de hasta 4.41×usando 4 hilos y solapamiento del 20 %. En el c´omputo total, se obtiene una aceleraci´on de hasta 7.87×. Recordemos que se podr´ıa conseguir una aceleraci´on mayor aumentando el n´umero de hilos que se usan en la detecci´on en el caso de que no fuera necesario registrar este tipo de escalas. Si comparamos los resultados entre ambas im´agenes, se puede ver que en la detecci´on de regiones se produce una aceleraci´on mayor en Baraboo Hills, debido a su mayor tama˜no, que mejora la eficiencia de la paralelizaci´on en esa fase. Sin embargo, en la descripci´on pasa justo lo contrario. Esto se debe a que Santa Barbara Front, a pesar de ser la imagen m´as peque˜na de las utilizadas, se obtienen una gran cantidad de regiones (ver Tabla IV). Este desacoplamiento entre el tama˜no y el n´umero de regiones que se detectan es lo que produce que Santa Barbara Front sea de las im´agenes que m´as aceleraci´on experimentan, como se vio en la Figura 3, a pesar de ser la m´as peque˜na. V. Conclusiones En este art´ıculo se ha presentado una versi´on paralela del algoritmo HSI-MSER de registro de im´agenes hiperespectrales basado en caracter´ısticas. Este algoritmo adapta MSER y SIFT para su uso en im´agenes hiperespectrales, aprovechando la informaci´on espectral de forma eficiente. La implementaci´on propuesta paraleliza el algoritmo haciendo uso los distintos n´ucleos disponibles en la CPU as´ı como de la GPU, obteniendo una considerable reducci´on de los tiempos de ejecuci´on cuando la comparamos con la versi´on secuencial original. La implementaci´on propuesta ha sido evaluada con 338 Jornadas SARTECO 2023 Tabla IV: N´umero de parejas de regiones, porcentaje de parejas correctas, n´umero de regiones, n´umero de escalas registradas correctamente y tiempos de ejecuci´on de la versi´on secuencial y CUDA + OpenMP 4 sin solapamiento y con 20 % de solapamiento. Imagen Secuencial CUDA + OpenMP 4 CUDA + OpenMP 4 + 20 % Jasper Ridge Parejas 253 187 236 Regiones imagen 1 32538 33379 37879 Regiones imagen 2 32625 33387 38128 Porcentaje de parejas correctas 0.742 0.778 0.786 Tiempo (s) 43.95 5.22 5.72 N´umero de escalas registradas 15 6 10 Santa Barbara Front Parejas 618 531 630 Regiones 1 27601 28578 32843 Regiones 2 31688 32618 37379 Porcentaje de parejas correctas 0.727 0.732 0.756 Tiempo (s) 33.70 3.89 4.28 N´umero de escalas registradas 15 11 15 Santa Barbara Line Parejas 977 825 988 Regiones 1 24420 25335 28861 Regiones 2 25703 26812 30261 Porcentaje de parejas correctas 0.773 0.764 0.774 Tiempo (s) 33.78 5.17 5.64 N´umero de escalas registradas 17 12 15 Baraboo Hills Parejas 133 110 137 Regiones 1 24413 25317 28876 Regiones 2 22223 23077 26261 Porcentaje de parejas correctas 0.756 0.771 0.769 Tiempo (s) 32.83 5.00 5.38 N´umero de escalas registradas 8 7 7 Crown Point Parejas 693 572 704 Regiones 1 27937 28764 33281 Regiones 2 27412 28283 32607 Porcentaje de parejas correctas 0.219 0.24 0.24 Tiempo (s) 36.11 5.01 5.49 N´umero de escalas registradas 17 13 12 Tabla V: Tiempos de ejecuci´on en segundos y aceleraci´on de la versi´on CUDA + OpenMP 4 con un 20 % de solapamiento respecto a la versi´on secuencial para Baraboo Hills. Fase Secuencial (s) CUDA + OpenMP 4 + 20 % (s) Aceleraci´on Leer las im´agenes 0.92 0.92 — Copiar im´agenes a GPU —0.28 — Selecci´on de bandas 0.32 0.12 2.68× Detecci´on de regiones 8.29 1.88 4.41× Descripci´on de regiones 14.99 1.02 14.66× Emparejamiento 8.30 0.40 20.94× Registro 0.01 0.01 — Total 32.83 5.38 6.10× im´agenes accesibles p´ublicamente para las que hay datos disponibles en la literatura. Para la evaluaci´on se han tenido en cuenta tanto la aceleraci´on respecto a la implementaci´on secuencial original como la capacidad de registro. Se ha comprobado que, si bien al ejecutar la versi´on paralela del registro se producen Jornadas SARTECO 2023 339