Clasificación Tisular en GPU: aceleración y optimizaciones

Departamento de Arquitectura de Computadores E.T.S. Ingeniería Informática TESIS DOCTORAL Clasificación tisular en GPU: aceleración y optimizaciones Presentada por: Antonio Ruiz Sánchez Dirigida por: Manuel Ujaldón Martínez Málaga, 2015 AUTOR: Antonio Ruiz Sánchez http://orcid.org/0000-0001-5997-9792 EDITA: Publicaciones y Divulgación Científica. Universidad de Málaga Esta obra está bajo una licencia de Creative Commons Reconocimiento-NoComercial-SinObraDerivada 4.0 Internacional: http://creativecommons.org/licenses/by-nc-nd/4.0/legalcode Cualquier parte de esta obra se puede reproducir sin autorización pero con el reconocimiento y atribución de los autores. No se puede hacer uso comercial de la obra y no se puede alterar, transformar o hacer obras derivadas. Esta Tesis Doctoral está depositada en el Repositorio Institucional de la Universidad de Málaga (RIUMA): riuma.uma.es Dr. D. Manuel Ujaldón Martínez. Profesor Titular del Departamento de Arquitectura de Computadores de la Universidad de Málaga. CERTIFICA: Que la memoria titulada “Clasificación tisular en GPU: Aceleración y optimizaciones”, ha sido realizada por D. Antonio Ruiz Sánchez bajo mi dirección en el Departamento de Arquitectura de Computadores de la Universidad de Málaga y concluye la Tesis que presenta para optar al grado de Doctor por la Universidad de Málaga. Málaga, a 4 de junio de 2015 Dr. D. Manuel Ujaldón Martínez. Director de la tesis. . A Indira, y a todos quienes la aman . Agradecimientos Esta tesis es el resultado del equilibrio de un cúmulo de circunstancias y sentimientos tan enfrentados como irracionales. Durante todo este tiempo he vivido muchos momentos de ilusión, de trabajo intenso, de abandonos, de aprendizaje, de decepción, de segundas oportunidades, etc. Sin embargo, tengo que dar las gracias a cada uno de estos momentos, que junto con el apoyo incondicional de mis más allegados, me han hecho ser mejor y crecer en lo personal y profesional. En primer lugar, quisiera expresar mi enorme gratitud a mi director de tesis, Manuel Ujaldón, por confiar y creer en mí incondicionalmente, por ser el mejor ejemplo del trabajo y el sobresfuerzo, y por transmitirme ánimos e ilusión incluso cuando ya daba por hecho que esta tesis quedaría en un cajón y nunca llegaría a ver la luz. Aún estoy más agradecido por sentir tu confianza incluso después de tomar la decisión de desvincularme de la investigación a tiempo completo. Por todo ello, esta tesis te pertenece, y me gustaría que la sintieras como un éxito de tu dirección y dedicación, a pesar de todas las tempestades ajenas a tu voluntad que has encontrado en el camino. Además, me gustaría resaltar y agradecer la oportunidad que me brindaste para trabajar en el área de las GPUs y la Supercomputación, tengo que reconocer que no he encontrado ningún tema que me llegue a apasionar tanto después de más de seis años de carrera profesional paralela. Tu enseñanza es un legado para afrontar con éxito mi carrera profesional. Me gustaría agradecer la hospitalidad de todo el personal del Departamento de Arquitectura de Computadores, desde becarios hasta profesores, sin olvidarme de los técnicos y la secretaria. Especial reconocimiento merece el interés mostrado por mi trabajo y las sugerencias recibidas del profesor Nicolás Guil. Quisiera hacer extensiva mi gratitud al personal del departamento de Biomedical Informatics de The Ohio State University. A los profesores Umit Catalyurek, Kun Huang y Metin Gurcan, a los ya doctores Tim Hartley, Jun Kong, Lee Cooper, y en especial a Olcay Sertel, quien me acogió durante mi estancia de investigación como a un hermano y con quien mantengo la amistad desde la distancia. II Agradecimientos Siento la obligación y necesidad de agradecer con mención individualizada a cada uno de mis compañeros del despacho de becarios (también conocido como corral). La obligación, porque ellos tuvieron que aguantarme en el día a día; la necesidad, porque además hicieron parecer que era un placer compartir tiempo y lugar. A todos ellos, Antonio Muñoz, Francisco Jaime, Ricardo Quislant, Rosa Castillo, Adrián Tineo, Juan Lucena, Sergio Varona, Miguel Ángel Sánchez, Victoria Martín, Javier Ríos, Maximiliano García, Fernando Barranco y Alfredo Martínez. Y a cada uno de mis compañeros de trabajo que se han involucrado de una manera u otra, especialmente a Pedro Antonio por su interés y recomendaciones tras la lectura de esta tesis. Quisiera también agradecer a cada uno de mis amigos de toda la vida, por el simple hecho de seguir siéndolos como el primer día. En especial a Salva, Escarcena, Chacón, Aroro y Trivi, ellos siempre han estado disponibles para todo. Por último, pero no por ello menos importante, me gustaría agradecer y dedicar este trabajo a mis padres, mujer e hija. Todos ellos han hecho realidad el llegar hasta aquí, con el apoyo incondicional y el seguimiento desde el silencio para encontrarme con todas las facilidades del mundo. Padres, gracias por vuestro apoyo y esfuerzo por permitirme realizarme en mis estudios y mi vida. Agradezco los consejos sabios que en el momento exacto habéis sabido darme para continuar luchando pese a todas las adversidades, y cuando pareciera que no los valoraba. A mi mujer, Ana Eva, debo agradecer encarecidamente todo su apoyo, anímico y afectivo, pese al tiempo que esta tesis le ha robado. Gracias por darle vida a la persona más importante de nuestras vidas. Indira nos ha enseñado lo que es amar de verdad. Gracias hija por tus ganas de jugar y disfrutar cuando veías que no te dedicaba todo el tiempo que merecías porque lo dedicaba a algo que aún estás lejos de entender. A cada uno de los mencionados en este párrafo, os quiero de una manera muy especial, ahora que soy padre sé de primera mano todo lo que siempre me aman los míos. Resumen Desde hace una década, los procesadores gráficos o GPUs vienen ganando protagonismo en la computación de altas prestaciones, contribuyendo a la aceleración de miles de aplicaciones en multitud de áreas de la ciencia. Pero más que esta conquista, lo que ha hecho singular al movimiento GPGPU ha sido la vía para su consecución, ofreciendo tecnología popular, barata y notablemente arropada. Como resultado, la supercomputación está hoy al alcance de cualquier usuario y empresa, democratizando un sector hasta entonces circunscrito a unos pocos centros elitistas. El auge de las GPUs en los entornos de altas prestaciones ha generado un reto a la comunidad de desarrolladores software. Los programadores están habituados a pensar y programar de manera secuencial, y sólo una minoría se atrevía hace 10 años a adentrarse en este mundo. La programación paralela es una tarea compleja que exige otras habilidades y modelo de razonamiento, además de conocer nuevos conceptos hardware, algoritmos y herramientas de programación. Poco a poco, esta percepción ha ido cambiando gracias a la aportación de aquellos que, conscientes de la dificultad, quisieron aportar su granito de arena para facilitar esta transición. El trabajo de esta tesis recoge este espíritu. Planteamos nuevos diseños e implementaciones de algoritmos en el ámbito de la biocomputación para evaluar el rendimiento de las GPUs más destacadas durante la última década, desde equipos con una única GPU hasta supercomputadores de 32 GPUs. En cada uno de los problemas de biocomputación se han analizado todas las características relevantes de la GPU que permiten exprimir su gran potencial, para así presentar de una manera didáctica y rigurosa un estudio pormenorizado de los detalles y técnicas de programación más acordes a cada tipo de algoritmo. Cronológicamente, la aparición de la arquitectura de cálculo paralelo CUDA para GPUs es un hito de especial importancia en la programación de algoritmos de propósito general en GPUs. Nuestro trabajo comenzó en la era pre-CUDA con una aplicación de detección de círculos basada en la transformada de Hough y un algoritmo de detección del tumor neuroblastoma. Sus implementaciones explotan la GPU desde una XÍNDICE GENERAL 7.2.5. Análisis de color y texturas sobre procesadores gráficos . . . . 142 7.2.6. Análisis del rendimiento de la arquitectura Kepler mediante algoritmos dinámicos e irregulares . . . . . . . . . . . . . . . 143 7.3. Línea de Investigación Actual y Trabajo Futuro . . . . . . . . . . . . 143 Bibliografía 145 Índice de figuras 1.1. El supercomputador BALE. . . . . . . . . . . . . . . . . . . . . . . . 11 2.1. Detalle del cauce de segmentación gráfico. . . . . . . . . . . . . . . . 14 2.2. Diagrama de bloques de la arquitectura Nvidia G80. . . . . . . . . . . 20 2.3. La interfaz hardware de CUDA para la GPU. . . . . . . . . . . . . . 21 3.1. Muestra representativa del proceso de la transformada Hough. . . . . 32 3.2. Primitivas OpenGL disponibles para que el rasterizador de la GPU haga la interpolación de una lista de vértices. . . . . . . . . . . . . . 34 3.3. Descripción de la transformada de Hough en GPU. . . . . . . . . . . 35 3.4. Proceso de interpolación sobre un único punto del contorno de los votos generados por el rasterizador de la GPU. . . . . . . . . . . . . . 36 3.5. Proceso de acumulación de votos sobre una textura gráfica. . . . . . . 37 3.6. Unidades funcionales del cauce de segmentación gráfico empleadas para la implementación de la THC en GPU. . . . . . . . . . . . . . . 38 3.7. Espacio parametrizado de la transformada Hough como resultado del rasterizador de la GPU sobre una textura. . . . . . . . . . . . . . . . 39 3.8. Resultados de las estrategias implementadas para la Transformada Hough para una imagen de 20 círculos con un radio de 50 píxeles. . . 40 3.9. Reducción del área de error al aumentar el número de semillas s. . . . 42 3.10. Cálculo de error producido cuando la transformada Hough se lleva a cabo a través de la primitiva GL_LINES para el caso particular de s=8. 44 3.11. Planos de renderizado empleados en las distintas optimizaciones en GPU. .................................. 48 XI XII ÍNDICE DE FIGURAS 3.12. Entradas y salidas de datos de la THC para las imágenes representativas c20r100 y c100r100. . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1. Muestras de imágenes de neuroblastoma. . . . . . . . . . . . . . . . 64 4.2. Diagrama de flujo del algoritmo de clasificación del estroma. . . . . . 65 4.3. Cálculo de la matriz de co-ocurrencia para una imagen de tamaño 4x4 en la que se muestra las intensidades de cada píxel. . . . . . . . . . . 66 4.4. Operador LBP sobre una imagen de tamaño de 3x3. . . . . . . . . . . 67 4.5. Factores de aceleración entre la versión C++ y Matlab para diferentes tamañosdeimagen............................ 73 4.6. Factores de aceleración entre la versión GPU y C++ para diferentes tamañosdeimagen............................ 74 4.7. Factores de aceleración entre la versión GPU y Matlab para diferentes tamañosdeimagen............................ 74 4.8. Matrices locales de co-ocurrencia en CUDA. . . . . . . . . . . . . . 80 4.9. Operador LBP en CUDA. . . . . . . . . . . . . . . . . . . . . . . . . 81 4.10. Estructura del DataCutter para la aplicación de análisis de imágenes. . 83 4.11. Comparativa de todas las implementaciones del análisis de imágenes sobre un nodo cuando la imagen es pequeña. . . . . . . . . . . . . . . 84 4.12. Comparativa de tiempos de ejecución de las implementaciones en GPU y DataCutter sobre un único nodo para los tres tamaños de imagen. 85 4.13. Tiempos de ejecución de C++, Cg y CUDA para la implementación conDataCutter. ............................. 87 4.14. Resumen del factor de aceleración para las diferentes implementacionesenparalelo. ............................. 87 5.1. Registro rígido rápido usando características de alto nivel. . . . . . . 94 5.2. Proceso de emparejamiento de características. . . . . . . . . . . . . . 96 5.3. Áreas seleccionadas aleatoriamente correspondiente al emparejamiento de características entre dos imágenes. . . . . . . . . . . . . . . . . 98 5.4. Flujo de trabajo del algoritmo de registro que consta de dos fases: el registro rígido y no rígido. . . . . . . . . . . . . . . . . . . . . . . . 100 5.5. Correlación cruzada normalizada basada en la FFT. . . . . . . . . . . 102 5.6. Demostración visual del registro realizado sobre algunas muestras. . . 106 ÍNDICE DE FIGURAS XIII 5.7. Porcentaje de características procesadas por imagen en cada conjunto de imágenes de entrada. . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.8. Tiempos de ejecución para el algoritmo de registro en una CPU Opteron.108 5.9. Tiempos de ejecución en la GPU Quadro para el algoritmo de registro. 109 5.10. Comparativa entre los tiempos de ejecución de la CPU y la GPU en términos de factores de aceleración. . . . . . . . . . . . . . . . . . . 110 5.11. Escalabilidad de la GPU: factores de mejora cuando se habilita una segundaGPU...............................112 6.1. Pseudocódigo para el cálculo de los momentos de Zernike. . . . . . . 118 6.2. Tiempo de ejecución en Kepler para evaluar el rendimiento en función deltamañodebloque...........................125 6.3. Ganancia obtenida con el uso de streams para distintos tamaños de imagen cuando aumentamos el conjunto de momentos de Zernike. . . 128 6.4. Beneficio cuando el tamaño de imagen corresponde a un bloque CUDA para aislar la aceleración atribuida a Hyper-Q y Kernels Concurrentes. .................................129 6.5. Comparativa entre la variante recursiva frente a Hyper-Q en el método directo. .................................132 Índice de tablas 1.1. Características de la CPU y GPU del PC utilizado en 2007. . . . . . . 9 1.2. Características hardware de las CPUs y GPUs usadas en cada nodo del supercomputador BALE. . . . . . . . . . . . . . . . . . . . . . . 10 1.3. Características de las GPUs del servidor de computación Yuca. . . . . 12 2.1. Principales rasgos CUDA de las GPU basadas en las generaciones FermiyKepler. ............................. 24 3.1. Lista de las primitivas OpenGL para unir vértices. . . . . . . . . . . . 34 3.2. Número de votos para las diferentes implementaciones en GPU. . . . 38 3.3. Error por voto con la primitiva GL_LINES en función del número de semillas y el tamaño del radio. . . . . . . . . . . . . . . . . . . . . . 45 3.4. Número mínimo de semillas, s, para mantener el error máximo por debajo de la distancia de un píxel para diferentes tamaños de radio. . . 47 3.5. Conjunto de optimizaciones empleadas con el rasterizador de la GPU en términos de precisión y eficiencia. . . . . . . . . . . . . . . . . . . 47 3.6. Tiempos de ejecución para diferentes imágenes y pasos de ángulo bajo dos CPUs diferentes: Intel Core 2 Duo y Pentium 4. . . . . . . . . 52 3.7. Tiempos de ejecución de la Transformada Hough para círculos de tamaño 50 píxeles de radio bajo diferentes estrategias. . . . . . . . . . . 53 3.8. Conjunto de imágenes empleadas en las pruebas experimentales. . . . 54 3.9. Votos emitidos e interpolados en función de la implementación GPU y el tamaño de radio para la transformada Hough de Círculos. . . . . . 56 3.10. Tiempos de ejecución y factores de aceleración para todas las optimizaciones desarrolladas para la transformada de Hough. . . . . . . . . 57 XV XVI ÍNDICE DE TABLAS 3.11. Tiempos de ejecución para diferentes estrategias y carga de trabajo de laTHC. ................................. 58 4.1. Precisión de clasificación para diferentes tamaños de la matriz de coocurrencia. ............................... 68 4.2. Tiempos de ejecución para una imagen de tamaño 1020 x 916 píxeles en un PC de escritorio del 2007. . . . . . . . . . . . . . . . . . . . . 71 4.3. Factor de aceleración para diferentes matrices de co-ocurrencia sobre una imagen de tamaño 1020 x 916 píxeles bajo diferentes plataformas. 72 4.4. Influencia de la carga de trabajo en las ganancias de rendimiento para cada una de las tareas involucradas en el análisis del neuroblastoma. . 73 4.5. Comparativa de la escalabilidad del hardware desde el lado de la CPU ylaGPU. ................................ 73 4.6. Precisión de los valores de salida para diferentes plataformas hardware. 75 4.7. Factores de aceleración obtenidos para una matriz de co-ocurrencia de 4x4 sobre diferentes tamaños de imagen. . . . . . . . . . . . . . . 76 4.8. Principales optimizaciones CUDA en la aplicación de análisis de imagen. ................................... 78 4.9. Tiempos de ejecución de la aplicación de análisis de imagen para diferentes métodos de programación y plataformas hardware. . . . . . . 82 4.10. Propiedades de las imágenes usadas en los experimentos. . . . . . . . 84 5.1. Peso porcentual medio de cada una de las fases antes y después de migrar el algoritmo hacia la GPU. . . . . . . . . . . . . . . . . . . . 100 5.2. Diferentes tamaños de las ventanas de búsqueda y características en los algoritmos de registro. . . . . . . . . . . . . . . . . . . . . . . . . 103 5.3. Conjunto de imágenes de entrada empleadas como datos de entrada para el algoritmo de registro. . . . . . . . . . . . . . . . . . . . . . . 104 5.4. Tamaños de las ventanas de búsqueda y características empleadas en los algoritmos de registro. . . . . . . . . . . . . . . . . . . . . . . . . 106 5.5. Tiempos de ejecución y factores de aceleración para el algoritmo de registro..................................109 5.6. Número de mosaicos procesados y descartados en una configuración condosGPUs. .............................112 ÍNDICE DE TABLAS XVII 6.1. Tiempos de ejecución para procesar todas las repeticiones de un orden a través del método directo de los momentos de Zernike. . . . . . . . 123 6.2. Tiempos de ejecución y factores de aceleración logrados para las diferentes estrategias que explotan el paralelismo dinámico. . . . . . . . 127 6.3. Tiempos de ejecución para todas las repeticiones de un orden específico a través del método q-recursive en las GPU Fermi y Kepler. . . . 130 Parte I Introducción 1 8 Capítulo 1. Aplicaciones biomédicas sobre arquitecturas gráficas El Capítulo 2 presenta en detalle las GPUs, desde las primeras con más de 10 años hasta las más recientes que hemos utilizado, basadas en la generación Kepler de Nvidia. La segunda parte está relacionada con la programación gráfica más clásica basada en la renderización como si se tratara de videojuegos o animaciones gráficas. Tanto la arquitectura gráfica como las instrucciones son equivalentes a las utilizadas en las animaciones gráficas, de modo que se delega en la pericia de los desarrolladores la transformación del problema computacional en otro ligado al cauce de segmentación gráfico. El Capítulo 3 implementa sobre esta perspectiva clásica un algoritmo para la detección de formas circulares en imágenes. Durante el transcurso de este capítulo se llevan a cabo, desde un puesto de vista didáctico, numerosas optimizaciones gráficas de la transformada de Hough para la detección de círculos. Se utiliza para ello la arquitectura gráfica convencional, con elementos como los procesadores de vértices y píxeles, el rasterizador y las unidades de blending, que son evaluados desde el punto de vista del rendimiento y la precisión. La tercera parte resuelve diferentes problemas sobre plataformas gráficas basadas en CUDA (Compute Unified Device Architecture). Desarrollado por la empresa Nvidia, fabricante de más del 75% de las GPUs actuales, CUDA define tanto el modelo de programación como la arquitectura subyacente y el interfaz para la programación de aplicaciones de propósito general (API, Application Program Interface). Las prestaciones y facilidad de uso ofrecidas por CUDA han ido en aumento desde su aparición, y los capítulos de esta parte analizarán todas ellas desde la perspectiva de nuestras aplicaciones biomédicas. El Capítulo 4 emplea CUDA para desarrollar un algoritmo de detección del tumor neuroblastoma. Este capítulo servirá de introducción a una versión inicial de CUDA, puesto que se implementó durante la transición entre la programación con Cg (C for graphics) y la programación con CUDA para aplicaciones de propósito general. Además, aprovechando el momento del cambio, se evaluará una interesante comparativa entre las diferentes metodologías de acceso a la GPU, desde equipos con una sola GPU hasta supercomputadores dotados de 16 nodos con 2 GPUs cada uno. El Capítulo 5 ilustra la potencia de cálculo de las primitivas proporcionadas en la librería de CUDA. En concreto, se emplea el cálculo de la transformada rápida de Fourier (FFT) en una aplicación de registro de imágenes. El Capítulo 6 aprovecha las últimas novedades de la arquitectura gráfica Kepler para implementar el algoritmo de los momentos de Zernike, muy útil para el reconocimiento de patrones en imágenes biomédicas, entre otras muchas aplicaciones. Este tramo final aprovecha ya las últimas prestaciones aparecidas en las GPUs de tercera generación correspondientes al bienio 2013-2014. 1.3. Descripción de los sistemas empleados durante los experimentos 9 La cuarta parte presenta las principales conclusiones y aportaciones de la investigación de esta tesis. La Sección 7.1 describe las conclusiones generales obtenidas de cada uno de los trabajos de investigación realizados para cada Capítulo. La Sección 7.2 muestra las diferentes publicaciones fruto del trabajo de esta tesis. Finalmente, la Sección 7.3 describe brevemente las líneas futuras de investigación que nacen de esta tesis. 1.3. Descripción de los sistemas empleados durante los experimentos Para evaluar las implementaciones presentadas en los siguientes capítulos se han utilizado plataformas de muy variado coste y complejidad, todas ellas con una arquitectura híbrida CPU-GPU, y representativas de la época en la que se abordaron dichas implementaciones. A continuación se describen en orden cronológico cada una de ellas. Características del procesador CPU 2007 GPU 2007 Fabricante Intel Nvidia Modelo Core 2 Duo GeForce 8800 GTX Arquitectura E6400 Conroe G80 Frecuencia 2.13 GHz 575/1350 MHz Potencia bruta de procesamiento 10 GFLOPS 520 GFLOPs Tipo de memoria DDR2 GDDR3 Tamaño 4 GBytes 768 MBytes Ancho del bus 2 x 64 bits 384 bits Frecuencia 2 x 333 MHz 2 x 900 MHz Ancho de banda 10.8 GB/s 86.4 GB/s Tabla 1.1: Características de la CPU y GPU del PC utilizado en 2007. 1.3.1. PC típico de 2007 con una CPU y GPU La primera configuración y más simple trata de un PC del 2007 compuesto de una única CPU y GPU. En la Tabla 1.1 se detallan las especificaciones de sus procesadores. La GPU con arquitectura G80 es la primera que incorpora la tecnología CUDA y que, aunque mantiene la compatibilidad con los algoritmos de propósito general GPGPU diseñados con Cg, proporciona un entorno más amigable al desarrollador. En los capítulos siguientes se analiza una comparativa para esta GPU entre la programación con Cg y con CUDA para los mismos algoritmos. 10 Capítulo 1. Aplicaciones biomédicas sobre arquitecturas gráficas 1.3.2. Supercomputador BALE del Ohio Supercomputer Center La segunda configuración de la misma época pero más enfocada a supercomputación es el clúster BALE del centro estadounidense Ohio Supercomputer Center. BALE está compuesto por un total de 71 nodos Linux, pero para nuestro propósito aquí lo más interesante se encuentra en los 16 nodos de visualización incorporados, cada uno de ellos equipados con dos CPUs de doble núcleo AMD Opteron 2218 y dos GPUs Nvidia Quadro FX 5600. La red de interconexión entre los nodos es Infiniband. La Figura 1.1 ilustra la arquitectura de los nodos de visualización, y la Tabla 1.2 resume las especificaciones de los dos procesadores, tanto CPU como GPU. Cada nodo de BALE es una placa base con doble zócalo (dual-socket) dotada de sendos procesadores Opteron X2 2218 de doble núcleo a una frecuencia de 2.6 GHz. Cada núcleo en el sistema tiene un par de cachés L1 gemelas para datos e instrucciones, con capacidad para 64 KB y asociatividad en conjuntos de 2 vías. La caché L2 de 1 MB y también asociativa de dos vías no es compartida por los núcleos, pero sí mantiene la coherencia de caché. Cada zócalo proporciona su propio controlador de memoria DDR2 a 667 MHz de doble canal al igual que un enlace HyperTransport para acceder a la caché y memoria de los otros zócalos, ofreciendo un ancho de banda de 10.6 GB/s para un total de 21.3 GB/s para cada nodo de 8 GB de memoria principal. El pico de rendimiento para una aritmética de doble precisión es de 4.4 GFLOPS por núcleo, 8.8 GFLOPS por zócalo y 17.6 GFLOPS por nodo. En simple precisión, el rendimiento conseguido por nodo asciende a 35.2 GFLOPS, proporcionando un total para todo el conjunto de nodos de visualización de 563.2 GFLOPs. Características CPU de AMD GPU de Nvidia Placa base ASUS KFN32-D SLI Quadro FX 5600 Modelo de procesador Opteron X2 2218 G80 Velocidad del procesador 2.6 GHz 600/1350 MHz Número de zócalos 2 2 Número de núcleos 2 128 Rendimiento pico 2 x 8.8 GFLOPS 2 x 330 GFLOPS Tamaño memoria 8 GBytes 2 x 1.5 GBytes Ancho del bus 2 x 64 bits 2 x 384 bits Frecuencia memoria 667 MHz 1600 MHz Ancho de banda 2 x 10.8 GB/s 2 x 76.8 GB/s Tabla 1.2: Características hardware de las CPU y GPUs usadas en cada nodo del supercomputador BALE. Los GFLOPS se calculan para aritmética de punto flotante de simple precisión (32 bits). Respecto a las prestaciones gráficas de cada nodo, tiene dos tarjetas gráficas Nvidia Quadro FX 5600 basadas en la arquitectura G80. En un entorno gráfico, la ar- 1.3. Descripción de los sistemas empleados durante los experimentos 11 Figura 1.1: El supercomputador BALE. quitectura G80 puede verse como un cauce de segmentación gráfico de 4 fases para sombreadores (shaders), texturas, rasterizado y coloreado. Como arquitectura paralela, la G80 es un procesador SIMD (Single Instruction Multiple Data) compuesto de 128 núcleos, accesible a través de la interfaz proporcionada por CUDA. 12 Capítulo 1. Aplicaciones biomédicas sobre arquitecturas gráficas 1.3.3. Servidor de computación YUCA La última y más moderna configuración la aporta el servidor de computación Yuca del Departamento de Arquitectura de Computadores de la Universidad de Málaga (ver Tabla 1.3). Las dos GPUs empleadas son diferentes: una con arquitectura Fermi (vigente en el trienio 2010-2012) y otra con arquitectura Kepler (2013-2015). La GPU Fermi permite situar la referencia de los resultados experimentales en una arquitectura de segunda generación, de modo que se pueda cuantificar la rémora con respecto a la tercera generación de multiprocesadores CUDA SMX, incluso antes de aplicar las nuevas técnicas y características introducidas por Kepler. Procesador GPU Fermi (Nvidia) GPU Kepler (Nvidia) Unidades 1 1 Modelo comercial Tesla C2075 Tesla K20c Números de núcleos @ frecuencia 448 @ 1.15 GHz 2496 @ 0.71 GHz Threads activos por núcleo 48 64 Rendimiento pico 1.03 TFLOPS 3.52 TFLOPS Frecuencia de la memoria 2x 1566 MHz 2x 2600 MHz Ancho de banda del bus 384 320 Ancho de banda de memoria 148 GB/s 208 GB/s Tamaño de memoria y tipo 6 GB de GDDR5 5 GB de GDDR5 Bus a/desde CPU PCI-e x16 2.0 PCI-e x16 2.0 Tabla 1.3: Características de las GPUs del servidor de computación Yuca. 2Arquitecturas gráficas La popularidad alcanzada por las GPUs como procesadores gráficos programables se debe fundamentalmente a su bajo coste y al gran rendimiento alcanzado en muchas aplicaciones de propósito general. Sin embargo, la naturaleza de las arquitecturas gráficas es bien diferente respecto a la actual evolución de estos procesadores dentro de la computación de altas prestaciones. Esta introducción a la arquitectura de la GPU se desarrolla en orden cronológico para conocer su origen y metamorfosis, lo que nos permitirá entender mucho mejor el modelo de programación implementado y su adecuación a ciertos tipos de algoritmos. Hemos estructurado este capítulo de la siguiente forma: La Sección 2.1 presenta la arquitectura gráfica en sus orígenes y describe sus unidades funcionales. La Sección 2.2 muestra el cambio de arquitectura gráfica para habilitar las tecnologías y características específicas de una evolución hacia la computación de propósito general. La Sección 2.3 introduce las características del modelo de programación para las arquitecturas gráficas actuales. Finalmente, la Sección 2.4 recoge la evolución de las características implementadas para mejorar el rendimiento de arquitecturas gráficas cuando éstas se orientan a propósito general. 2.1. Cauce de segmentación gráfico La naturaleza de la propia arquitectura gráfica justifica que el diagrama de bloques en sus orígenes mantuviera una composición y distribución de componentes común independientemente del fabricante. Los elementos que constituían una tarjeta gráfica de mediados de los noventa siguen ahí, tan sólo levemente retocados, y son los siguientes: 13 14 Capítulo 2. Arquitecturas gráficas El procesador gráfico (GPU), encargado de mover los vértices iniciales y agruparlos en polígonos, para posteriormente rasterizarlos y transformarlos en píxeles aplicando texturas y colores, que finalmente se mezclarán con los objetos ya existentes en la escena. En el argot gráfico, este proceso recibe el nombre de renderización, y sigue vigente en esencia, aunque las transformaciones de vértices y píxeles son programables mediante sombreadores (shaders) (2001), cuya polivalencia y unificación originaron CUDA (2006). La Figura 2.1 ilustra el proceso en su conjunto. Figura 2.1: Detalle del cauce de segmentación gráfico. La memoria de vídeo (VRAM), que aloja todos los datos necesarios para llevar a cabo la renderización, desde los vértices y sus atributos al comienzo del cauce de segmentación, hasta su presentación final en el frame buffer, que no es más que una representación interna de lo que vemos en pantalla. Además de estos dos elementos, se alojan en memoria de vídeo los mapas de texturas y el zbuffer1obuffer de profundidad. Los procesadores de señal, luz y color (RAMDAC), cuya principal misión es convertir la información digital almacenada en el frame buffer en señales aptas para los dispositivos externos de vídeo, analógica para VGA o digital en fechas más recientes para DVI o HDMI. Los conectores externos, que permiten obtener el resultado de la información 1El buffer de profundidad que permite distinguir la región visible y oculta de cada objeto según la tercera dimensión espacial (plano Z). 2.1. Cauce de segmentación gráfico 15 procesada en diferentes fuentes o formatos. Entre ellos podemos citar como más populares DVI, VGA, HDMI, Super-VHS, HD-TV y vídeo compuesto. La BIOS de vídeo. Es el firmware de configuración de la tarjeta gráfica. Esta función ha ido delegándose hacia los drivers de la tarjeta gráfica, existiendo utilidades que permiten manipular la configuración desde del propio sistema operativo (SO). Los buses de comunicación, siendo el bus de memoria de vídeo y el bus que conecta el zócalo gráfico con la CPU (AGP o PCI-Express) los dos más relevantes. El primero comunica la GPU con su memoria de vídeo, y el segundo con el procesador central o CPU, normalmente a través del puente norte, el chip de la placa base que implementa sus controladores más importantes. Los disipadores encargados de refrigerar todo el conjunto y mantener una temperatura de trabajo adecuada. Suelen ser los más aparatosos y los que ocultan el resto de elementos ya mencionados. 2.1.1. Unidades Funcionales de la GPU El elemento principal y diferenciador, además de ser el más complejo en cuanto a su arquitectura, es la GPU. De hecho, su protagonismo lleva al lenguaje de la calle a decir que una tarjeta gráfica es de Nvidia o ATI, cuando en realidad sólo lo es la GPU. Nvidia no comercializa GPUs, aunque sí ensambla unas pocas para donarlas a sus GPU Education Centers y GPU Research Centers. El Departamento de Arquitectura de Computadores de la Universidad de Málaga posee ambas distinciones desde hace unos años, lo que le ha permitido beneficiarse de multitud de donaciones, y así las tarjetas gráficas que hemos utilizado en esta tesis proceden directamente de la firma estadounidense. La complejidad de la GPU ha crecido sobremanera, alcanzando e incluso, en algunos casos, superando a los procesadores de propósito general. Sus elementos más importantes son los procesadores de vértices y píxeles, responsables en buena parte de su versatilidad. Cada unidad funcional trata la información recibida y la devuelve procesada a la siguiente etapa del cauce de segmentación gráfico. A continuación desglosamos de forma más pormenorizada cada una de ellas. 16 Capítulo 2. Arquitecturas gráficas Vértices Los vértices son la base de la representación de la geometría de los objetos y contienen múltiple información en forma de atributos: Coordenadas tridimensionales. Cada vértice representa su posición por un vector de cuatro componentes. Las dos primeras (x, y)representan la posición en el espacio de dos dimensiones, la tercera zla profundidad o la posición que ocupa en el z-buffer, y la última, w, es un factor que permite normalizar el vector (x, y, z, w) = (x/w, y/w, z/w, 1). Color. Los vértices pueden presentar dos tipos de colores, el real y el reflectante para los efectos de luz. Ambos tipos son definidos por otras cuatro componentes: rojo, verde, azul y alfa. Sigue el modelo de color RGB con un canal adicional, alfa, para controlar el grado de transparencia u opacidad de cada píxel. Primitivas Una primitiva es una forma geométrica simple, desde un punto hasta un polígono. El tipo de primitiva más usado en computación es el triángulo por su versatilidad y rendimiento. El uso de vértices permite formar diferentes formas de primitivas. En la Sección 3.2 hay un ejemplo práctico del uso de diferentes primitivas y su descripción detallada. Agrupación de vértices Esta tarea, además de asociar los vértices con las primitivas seleccionadas, permite generar más vértices con el fin de aumentar el nivel de detalle de la superficie. Para ello, la GPU necesita las posiciones de los vértices y sus vectores normales. Existen varios métodos para realizar esta tarea, el más simple consiste en aumentar el número de vértices hasta llegar a la resolución deseada, con la desventaja de tener los objetos más distantes con más detalle de lo apreciable y los cercanos pixelados. El método ideal consiste en aumentar el número de vértices conforme disminuye la distancia al punto de vista (z-buffer). 2.1. Cauce de segmentación gráfico 17 Procesado de vértices El procesado de vértices es uno de los elementos más importantes de la GPU. Este procesador recibe los vértices y los trata en una de las dos modalidades permitidas: fija y programable (shader de vértices). La funcionalidad fija es la única permitida desde los orígenes de las tarjetas gráficas, cuya misión principal es la de efectuar transformaciones basadas en traslaciones, desplazamientos, giros y cambios de escala sobre el conjunto de los vértices. La funcionalidad programable permite procesar los vértices con programas específicos, denominados shaders, que permiten expresar ricos efectos en un lenguaje definido como tal. Este lenguaje ha evolucionado desde HLSL (High Level Shading Lenguage) hasta Cg (C for graphics), que puede considerarse una antesala de CUDA con un nivel de abstracción superior. Clippling,culling y rasterización Los vértices tratados necesitan ser acotados con diversas técnicas para representarlos en un plano bidimensional con unas dimensiones predeterminadas. La operación clipping descarta los polígonos que quedan fuera del ángulo de visión definido por el punto de vista y el plano de la imagen. Define por tanto los límites del espacio bidimensional a representar. La operación culling descarta los vértices de cada objeto que quedan ocultos tras su cara visible. Ambas operaciones persiguen filtrar el número de vértices que tienen relevancia para componer la imagen, aliviando la carga computacional de las tareas que les suceden. La etapa de rasterización convierte cada punto, línea o polígono 3D en una matriz 2D de puntos donde se guarda la información acerca del color y la profundidad para cada uno de los puntos que lo conforman. Para ello entra en juego la resolución de pantalla, ubicándose los vértices en ella mediante un proceso de interpolación. Como resultado, la lista de vértices de entrada se ha convertido en una matriz cuyas dimensiones corresponden a las del frame buffer. Procesado de píxeles o fragmentos El procesador de píxeles o fragmentos transforma los píxeles atendiendo a la información que les acompañan en forma de atributos. Un buen ejemplo de estos atributos son las coordenadas para acceder al mapa de texturas que se quiere utilizar. Al igual que el procesador de vértices, éste también dispone de una funcionalidad fija 24 Capítulo 2. Arquitecturas gráficas Generación de GPU Fermi Kepler Modelo hardware GF100 GK110 Hilos por warp 32 32 Máximo número de warps por multiprocesador 48 64 Bloques activos por multiprocesador 8 16 Máximo tamaño de bloque (en hilos) 1024 1024 Máximo número de hilos por multiprocesador 1536 2048 Máximo número de registros por hilo 63 255 Dimensión máxima de la malla de hilos 216 −1 232 −1 Paralelismo dinámico No Sí Hyper-Q No Sí Tabla 2.1: Principales rasgos CUDA de las GPU basadas en las generaciones Fermi y Kepler. cualidades, el paralelismo dinámico y la planificación Hyper-Q ofrecen gran potencial como alternativa en computación irregular, aunque eso sí, de forma más exigente para el programador. 2.4.1. Paralelismo dinámico En un sistema híbrido CPU-GPU, la ejecución eficiente de aplicaciones con elevado grado de paralelismo depende en gran medida de la versatilidad de los mecanismos que permitan distribuir el trabajo entre ambos aprovechando las mejores cualidades de cada uno de ellos. Hasta 2013, CUDA postulaba a la GPU como un coprocesador esclavo de la CPU que recibía sus encargos y trataba de acelerarlos, pero con una autonomía bastante limitada. El paralelismo dinámico permite a la GPU lanzar sus propios kernels, crear los eventos e hilos necesarios para controlar las dependencias, sincronizar los resultados y controlar la planificación de tareas, todo ello sin intervención alguna de la CPU, que ahora puede dedicarse de forma más eficiente a sus propias tareas. Por su parte, la GPU aporta un procesamiento más directo de bucles anidados y algoritmos recursivos, y en general, una computación más natural de código dinámico y estructuras de datos irregulares. Por ejemplo, ahora es posible determinar en tiempo de ejecución el número de hilos encargados de procesar los nuevos kernels creados, pudiendo establecerse una configuración inicial con un paralelismo más conservador que evite cálculos innecesarios en zonas livianas, para aumentar gradualmente el paralelismo en aquellas zonas que vayan requiriendo una computación más exigente. 2.4. Evolución de las características y funcionalidades en la arquitectura CUDA 25 2.4.2. Hyper-Q En el modelo CUDA, la CPU lanza los kernels sobre la GPU de forma secuencial, esto es, toda la GPU se dedica a procesar el primer kernel y hasta que éste no concluye no se inicia el segundo, y así sucesivamente. En los casos en que los kernels sean independientes y no quiera articularse esta barrera implícita de sincronización, puede utilizarse el concepto de streams para agrupar los kernels dependientes en un mismo stream y separar los kernels independientes en distintos streams. La búsqueda de un planificador óptimo que administre la carga de trabajo en GPU cuando ésta procede de diferentes streams es uno de los retos más difíciles de su arquitectura. Fermi permite una ejecución concurrente de hasta 16 streams, pero la existencia de una sola cola de trabajo obliga a multiplexar los streams, y por tanto, a su serialización. Aunque esta dependencia puede ser aliviada en una primera fase reordenando los kernels de cada stream, la tarea empieza a ser complicada y el rendimiento disminuye conforme la complejidad de los programas aumenta. Hyper-Q habilita hasta 32 colas de trabajo entre el host y el distribuidor de trabajo de CUDA en GPU, dotando de gran flexibilidad al conjunto para lograr grandes mejoras sin modificar la implementación. Ahora, cada stream se gestiona desde su propia cola hardware de trabajo, sin interferir las dependencias con otros streams, que pueden proceder del mismo programa CUDA u otros ubicados en diferentes procesos MPI o hilos POSIX (más conocidos como p-threads). De esta forma, la concurrencia es natural y no requiere preprocesamiento. Además, este mecanismo resulta más potente a medida que aumentamos el número de núcleos de cada multiprocesador, erigiéndose en uno de los pilares para la escalabilidad de las futuras generaciones de GPUs. Parte II GPGPU clásica: Programando gráficos 27 3Detección de Patrones La biología celular dispone de multitud de aplicaciones resueltas a través de un análisis detallado de imágenes previamente procesadas con alguna tinción química que resalte la estructura de los tejidos. Un ejemplo es la regeneración esquelética objeto de este capítulo, donde se pretende descubrir el microambiente que mejor induce la expresión de proteínas de matriz esquelética, responsables de diferenciar entre cartílago y hueso inducido a partir de células madres mesenquimales pluripotentes inyectadas en la región objeto de una lesión ósea. La evaluación del fenotipo de las imágenes tratadas con células para controlar los cambios producidos en el organismo se lleva a cabo con el recuento de células en cada tipo de tejido: cartílago, hueso y fibra. Este proceso pretende ser semi-automático para que el experto en la materia pueda apoyar su criterio en unos parámetros computables. De otra manera, el procedimiento se torna tedioso, sensible al observador y tremendamente costoso. Entre las técnicas de reconocimiento de patrones, que pueden ser clasificadas en dos familias dependiendo de si se emplea un algoritmo local o global, se ha seleccionado la transformada de Hough por ser una técnica de reconocida robustez, incluso con la presencia de solapamiento de formas y ruido. En este capítulo, se propone la implementación en GPU de la transformada de Hough para la detección de círculos junto con una serie de optimizaciones que permiten mejorar el rendimiento y maximizar el uso de las unidades funcionales del cauce de segmentación clásico de las arquitecturas gráficas. Esta implementación pone de manifiesto el potencial que ofrecen los procesadores de vértices, en desuso en una programación de propósito general en GPU, para problemas con unas determinadas características. Con el beneficio adicional de situarse en la fase inicial del cauce de segmentación de la GPU, permitiendo en una misma 29 30 Capítulo 3. Detección de Patrones fase de renderizado aprovechar otras unidades funcionales como los procesadores de píxeles y las unidades de blending (mezclado). El contenido de este capítulo se estructura de la siguiente forma. La Sección 3.1 describe el algoritmo completo para la detección de círculos así como el preprocesado necesario. La Sección 3.2 presenta el rasterizador de la GPU, que va a jugar un papel fundamental en las optimizaciones de la implementación en GPU de la transformada de Hough para detectar círculos (THC). Esta transformada se describe en la Sección 3.3, mientras que la Sección 3.4 analiza el impacto de las características de la textura usada como estructura de datos de entrada y salida. La Sección 3.5 estudia los parámetros que afectan a las diferentes implementaciones, mientras que la Sección 3.6 evalúa la precisión. Finalmente, las Secciones 3.7 y 3.8 presentan una serie de optimizaciones para GPU y CPU, respectivamente, cuyos resultados se exponen en la Sección 3.9. Finalmente, la Sección 3.10 resume las principales contribuciones de este capítulo. 3.1. El algoritmo La transformada Hough [5] es un método ampliamente utilizado para detectar círculos y contornos de diversas formas en el campo del procesamiento de imágenes. En la transformada Hough se parte de un objeto o figura objetivo y una imagen donde ha de realizarse la búsqueda del objeto. Posteriormente, se analizan todas las coincidencias del objeto en la imagen, que son caracterizadas como votos en un espacio de parámetros que representa la probabilidad de que el objeto exista en la imagen en una localización particular. Una tarea tan sencilla para el ojo humano resulta dificultosa y compleja por no disponer a priori de un patrón predefinido para todos los contornos en una imagen. La complejidad y el uso de memoria [24] se disparan a tal nivel que convierte dicho problema en fiel candidato para estudiar su paralelización. La automatización de dicha tarea se ha abordado desde diversas arquitecturas paralelas, como los sistemas multiprocesadores con memoria distribuida [26, 43, 78], computadores piramidales [3], máquinas SIMD [45], hardware de propósito específico [12] y arquitecturas reconfigurables [59, 61]. De carácter más reciente, la GPU ha sido igualmente empleada por su gran aceptación para la resolución de problemas exigentes computacionalmente [58]. Aunque el reconocimiento de círculos recae en la transformada Hough, previamente se aplica un preprocesado para aislar los puntos de interés o contornos en la imagen, eliminando así el ruido de las imágenes de entrada. 3.1. El algoritmo 31 3.1.1. Detección de contornos: Operador Canny Los métodos de detección de contornos utilizan los operadores gradiente como herramienta básica de procesado [65]. Entre todas las posibilidades, el operador Canny [14] está consolidado como una gran alternativa por sus excelentes resultados en imágenes con escala de grises; para imágenes en color, se requiere una conversión inicial aplicando el operador luminancia. El operador Canny es un procedimiento con varias fases que se detallan a continuación: 1. Primera convolución. Los contornos de la imagen son suavizados al aplicar el operador gaussiano que controla la cantidad de detalle en los contornos a la vez que elimina el ruido. 2. Segunda convolución. El operador Sobel aplicado en una imagen bidimensional y basado en la primera derivada, detecta los contornos gracias al gradiente de la intensidad de cada píxel. 3. Eliminación de los no-máximos. Los contornos detectados cuya intensidad no coincide con los máximos locales del gradiente se eliminan para disponer de contornos limpios. 4. Histéresis de doble umbral. Los valores de intensidad considerados en la imagen final son acotados a través de dos umbrales. Con dos umbrales, T1 >T2, los contornos aceptados serán aquellos que dispongan de algún píxel con un valor de intensidad por encima del umbral T2 y los demás valores por encima de T1. 3.1.2. Detección de patrones: La transformada Hough Una vez aplicado el filtro de Canny, los puntos o píxeles que describen los contornos forman los datos de entrada para la transformada Hough. Cada punto genera la figura a detectar como votos en un espacio de parámetros en el cual los máximos de la acumulación de todos los votos desvelan las coordenadas de las figuras detectadas (ver Figura 3.1). En el caso particular de la transformada Hough para la detección de círculos [43], en el que se centra este capítulo, cada punto genera una superficie tridimensional que describe la siguiente ecuación: (x−a)2+ (y−b)2=r2(3.1) Donde (x, y)son las coordenadas de cada punto, y (a, b)yrrepresentan el centro del círculo y el radio, respectivamente. Basando la detección en círculos con un radio 32 Capítulo 3. Detección de Patrones (a) Imagen de muestra (b) Detección de contorno (c) Salida con GL_POINTS y δt = 0,01 (d) Salida con GL_POINTS y δt = 0,1 Figura 3.1: (a) Imagen de muestra. (b) Detección de contornos y entrada a la transformada Hough después de aplicar el operador Luminancia y el filtro Canny. (c-d) Salida de la GPU donde el frame buffer acumula todos los votos del espacio parametrizado bidimensional. predefinido r, el proceso global empieza con una transformación entre un espacio de coordenadas (x, y)y el espacio de parámetros (a, b)discretizado, donde cada punto se convierte en el centro de una circunferencia de radio rque incrementa el número de votos por cada celda que pasa en el espacio de parámetros (a, b). Desglosando la ecuación 3.1 en cada eje, la conversión entre ambos espacios vendría definida por el siguiente par de ecuaciones: a=x−r∗cos(t) b=y−r∗sin(t)}t∈[0,2∗π](3.2) Con la descripción dada, es necesario en una fase inicial decidir la discretización para el ángulo ty el radio rde acuerdo a las características de la imagen. Conforme aumente la granularidad del ángulo, la carga computacional y la precisión aumentan, así como el tiempo de ejecución y la calidad de la detección de los círculos que son detectados en el espacio de parámetros como aquellos puntos votados con un valor superior a un umbral determinado. La disminución del paso de ángulo, δt, no siempre implica una mejor detección de círculos. Conforme disminuye la distancia entre votos, el peso porcentual de cada voto disminuye al igual que se distribuyen de una forma difuminada alrededor del centro real del círculo, dificultando la localización de un único punto como máximo local (ver Figura 3.1c con un δt superior a la Figura 3.1d). Existen métodos basados en la información del gradiente para la THC con unas exigencias más relajadas para la carga computacional y la memoria [21] pero que desafortunadamente sufren de imprecisión conforme el ruido es mayor [20]. El algoritmo 3.2. Impacto del rasterizador de la GPU 33 clásico es menos sensible al ruido a costa de una mayor carga computacional, y el objetivo de este trabajo es desarrollar una implementación mucho más eficiente en GPU que permita aprovechar la robustez del algoritmo en tiempos de ejecución razonables. 3.2. Impacto del rasterizador de la GPU Los datos de entrada para el cauce de segmentación gráfico consisten en un conjunto de comandos, texturas y vértices que son mandados desde la CPU a través del bus PCI-express. Haciendo uso de la librería gráfica OpenGL, los comandos inicializan y modifican estados, configuran renderizados y referencian vértices y píxeles; los vértices y los estados fluyen en la GPU a través de una secuencia de fases del cauce de segmentación (ver Figura 2.1): 1. El procesador de vértices permite la programación de cada vértice a través de un programa dedicado, denominado shader osombreador de vértices. 2. Los vértices son agrupados en primitivas que pueden definir diferentes formas, bien puntos, líneas o triángulos. Los bloques de Cull,Clip ySetup optimizan las operaciones de cada primitiva, eliminando caras no visibles, recortando la parte de la imagen no mostrada y desarrollando la configuración de las ecuaciones para los planos y contornos, respectivamente. 3. El rasterizador expande los datos hacia unos valores interpolados para generar todos los píxeles que cubren cada primitiva, dando un conjunto de fragmentos como resultado [49]. 4. Los píxeles candidatos (fragmentos) son aceptados por el procesador de fragmentos donde se aplica otro shader para realizar las transformaciones oportunas a cada píxel. 5. Los fragmentos resultantes son utilizados para escribir los resultados finales en el frame buffer (un array bidimensional de píxeles presentados en pantalla) o en una textura 2D como alternativa de almacenamiento. La lista de primitivas que puede emplear el rasterizador de la GPU para interpolar líneas y completar áreas está resumida en la Tabla 3.1 para el caso de OpenGL, con su funcionalidad mostrada en la Figura 3.2. Para el propósito de implementar la transformada Hough en GPU, las primitivas con relleno de área son descartadas y, por tanto, la selección se limita entre GL_POINTS y GL_LINES. GL_POINTS genera 40 Capítulo 3. Detección de Patrones (a) Espacio de parámetros ideal. GL_POINTS con δt =0.02. (b) Reducción de la carga computacional de la GPU. GL_POINTS con δt =0.1. (c) Uso del rasterizador de la GPU. GL_LINES con s= 16. Figura 3.8: Salidas de las diferentes estrategias implementadas para la THC usando una imagen de 1024x1024 píxeles que contiene 20 círculos con un radio de 50 píxeles. Las imágenes de la izquierda reflejan el potencial del frame-buffer a la hora de mostrar una salida visual que pueden orientar a los usuarios en las aplicaciones de reconocimiento semi-automático de círculos. Las gráficas de la parte derecha representan el espacio parametrizado alrededor de las coordenadas del centro de un círculo detectado; además una matriz de 3 x 3 muestra el porcentaje de los votos acumulados en los ocho píxeles vecinos al centro del círculo (el máximo valor de votado es de 100%). 3.4. Discretización y resolución de la textura 41 3.4.2. Déficit de votos Dada la discretización en GPU orientada hacia una textura bidimensional formada por píxeles, cuando un círculo tiene un radio de 10 píxeles, su perímetro, 2πr, estará formado por 2π10 = 62,83 píxeles. Este valor será redondeado a 63 o truncado hacia 62 píxeles, pero además el perímetro en ningún caso medirá el valor exacto proporcionado por la fórmula porque la distancia entre dos píxeles consecutivos difiere dependiendo de la alineación (horizontal o vertical) y, en caso de que estén situados en zonas del contorno donde forman una diagonal, aún estarían más lejanos y la contribución al perímetro sería mayor. Este déficit de votos encuentra su caso peor para la búsqueda de grandes círculos en resoluciones de textura más pobres. En la práctica, como una resolución mayor ofrece más píxeles a cada círculo, aumentar el tamaño del círculo o de la resolución de la textura ofrecen resultados similares en este sentido. Por ejemplo, el radio de la Figura 3.7b contiene 16 píxeles y el círculo completo está compuesto por un total de 92, cuando el valor teórico es de 2π16 = 100,53 píxeles (déficit del 8%); por otro lado, el radio de la Figura 3.7f es de 32 píxeles y el círculo está compuesto por 184, que difiere del valor teórico 2π32 = 201,06 píxeles en un 8%. El déficit esperado en el proceso de discretización puede ser moderado con la activación del modo antialiasing que suaviza los contornos con el uso del rasterizador. Para situaciones que requieran texturas de alta resolución, pueden aplicarse estrategias multiresolución para acelerar los tiempos de ejecución. La ejecución se daría desde un nivel inicial de baja resolución hasta una resolución final fijada, de modo que en cada resolución se retiene información valiosa para lograr la detección. En este sentido, la implementación que planteamos para la transformada Hough en GPU es perfectamente válida para poder elegir entre relajar la carga computacional o lograr un mayor grado de precisión. 3.4.3. Efectos colaterales La resolución de la textura afecta a otros aspectos de los métodos que usan las primitivas GL_POINTS y GL_LINES para calcular los votos en GPU: 1. Una baja resolución produce una deficiencia de votos. Se puede aprovechar la discretización de la textura en píxeles para seleccionar un incremento de ángulo, δt, que sólo genere votos en píxeles diferentes de acuerdo a la Ecuación 3.2, para así relajar la carga computacional sin afectar a la fidelidad de los resultados. 2. Una gran discretización para el ángulo puede provocar la pérdida de votos cuando se usa GL_POINTS. Para una combinación sensata entre buena precisión y 42 Capítulo 3. Detección de Patrones baja carga computacional, se puede querer mantener la resolución de la textura mientras se incrementa la discretización del ángulo para producir puntos más dispersos en la textura bajo la primitiva GL_POINTS (ver Figura 3.7g y 3.7h para incrementos de δt de 0.1 y 0.5). La Sección 3.5 comprueba que la discretización puede ser más agresiva en círculos pequeños. 3. Un bajo número de semillas para GL_LINES repercute en la precisión de los votos. La Figura 3.7c. muestra este efecto para una serie de semillas s = 32, 16, 8, 4; las líneas generadas por el proceso de interpolación en GL_LINES recuperan los votos situados entre las semillas calculadas. Dicha recuperación es menos precisa conforme disminuye el número de semillas. El papel de s es analizado en la Sección 3.5 desde el punto de vista de la precisión, y la Figura 3.9 muestra, en color negro, el área de error para la serie s = 4, 8, 16, 32. Las resoluciones más pequeñas afectan de igual manera tanto a los círculos procesados como a los interpolados. Figura 3.9: Reducción del área de error al aumentar sen un factor 4x desde su valor inicial hasta un valor de 32. El área de error se calcula dividiendo el área circular, πr2, por el área del polígono inscrito, sr2sin(π s) cos(π s). 3.5. Parámetros clave Desde una perspectiva cuantitativa, los parámetros determinantes para las implementaciones en GPU usando las primitivas GL_POINTS y GL_LINES son el incremento del ángulo δt y el número de vértices por círculo, respectivamente. En la Sección 3.6 se hace un estudio desde un punto de vista cualitativo. 3.5. Parámetros clave 43 3.5.1. El papel de δt Para la selección de un δt óptimo para la transformada Hough a través de la primitiva GL_POINTS, tiene que existir una correspondencia unívoca entre los votos generados por la Ecuación 3.2 y la posición de los puntos del contorno. Con una resolución suficiente para garantizar esta condición, el número de votos bajo una cierta discretización δt es de 2π δt , y el número de píxeles en el espacio de discretización para un círculo de radio rdebería ser de 2πr. El valor calculado bajo estas premisas es de δt =1 r, que revela una relación inversamente proporcional entre el tamaño de los círculos a detectar y la discretización del ángulo. De esta forma, la resolución del ángulo y la carga computacional se intensifican con el tamaño del círculo a detectar por la transformada Hough. El lado izquierdo de la Tabla 3.2 muestra algunos valores de δt con el número de votos procesadores para círculos de radio r = 10, 20, 50, 100. 3.5.2. El papel de s La parte derecha de la Tabla 3.2 muestra la influencia positiva del incremento de s para el beneficio del recuento de votos cuando los círculos en el espacio de parámetros son creados a través de interpolación con la primitiva GL_LINES. Aunque el número exacto de votos de la circunferencia debe ser el tamaño de la circunferencia 2πr, el número real para los círculos interpolados corresponde con el sumatorio del tamaño de los lados de los spolígonos regulares que describen la circunferencia, Perímetro = (Número de lados)·(Tamaño de cada lado) = s·2y(3.3) La Figura 3.10 representa el caso particular de s= 8, donde el tamaño de cada lado, 2y, puede calcularse mediante la siguiente expresión: r2=y2+x2⇒r2=y2+ (r·cos(π s))2⇒ y=√r2−(r·cos(π s))2⇒y=r·√1−cos2(π s) (3.4) Con scomo el número de lados y 2ycomo tamaño de cada lado, el perímetro del polígono proporciona la siguiente cantidad de votos: #votos =Perímetro polígono =s·2y= 2 ·r·s·√1−cos2(π s)(3.5) 44 Capítulo 3. Detección de Patrones Esta expresión refleja el déficit en el número de votos. Un aumento del déficit al incrementar el tamaño de la circunferencia puede ser compensado con el incremento del número de semillas, s(ver Figura 3.9). Al igual que en el análisis de δt, la resolución seleccionada deber ser suficiente para minimizar el déficit debido a las diferentes longitudes por la localización entre píxeles vecinos (vertical, horizontal y diagonal). Figura 3.10: Cálculo de error producido cuando la transformada Hough se lleva a cabo a través de la primitiva GL_LINES para el caso particular de s=8. 3.6. Precisión de la interpolación Como GL_POINTS garantiza una correcta localización de cada punto en el espacio de parámetros, el número máximo de votos que una textura puede alojar supone la restricción más importante para seleccionar el δt óptimo. Sin embargo, aunque s repercute en el número de votos usando GL_LINES, el proceso de interpolación proporciona una cantidad que difiere del valor teórico y que se acentúa con un número menor de semillas. El resto de esta sección analiza la precisión del rasterizador en el proceso de interpolación de votos. Por otro lado, la resolución de la textura del espacio de parámetros supone otra componente para la precisión como en cualquier proceso de discretización del dominio espacial. Sin embargo, para el estudio de la precisión se asume que la resolución de la textura es óptima, dicho valor se puede cuantificar como el error entre la diferencia de las áreas descritas por la población ideal de votos y los votos reales (ver Tabla 3.3). 3.6. Precisión de la interpolación 45 Radio - Error medio s = 4 s = 8 s = 16 s = 32 r = 10 píxeles 2.018 0.511 0.096 0.032 r = 20 píxeles 4.036 1.022 0.192 0.064 r = 50 píxeles 10.090 2.557 0.480 0.160 r = 100 píxeles 20.180 5.114 0.960 0.321 Error medio por voto 20.18 % 5.11 % 0.96 % 0.32 % como porcentaje de r Tabla 3.3: Tamaño medio del error por voto (en píxeles) para un radio rde un posible centro del círculo cuando la votación en el espacio de parámetros se realiza con GL_LINES con un número concreto de semillas, s. 3.6.1. Estimación del error El error total para la fracción perteneciente a un polígono de los sformados para generar el círculo corresponde a la diferencia entre el área de dicha porción del círculo y la recogida por el polígono formado (ver Figura 3.10), Error total =Error área = (área sector circular)−(área triángulo) =π·r2 s−r2·sin(π s)·cos(π s) = r2·(π s−sin(π s)·cos(π s)) (3.6) El error por voto es nulo justo en las semillas (como D en Figura 3.10), y alcanza su máximo valor en el punto equidistante entre dos semillas contiguas (como E en la Figura 3.10); dicho valor máximo se puede cuantificar como la diferencia entre el radio del círculo y la apotema de los polígonos generados usando ssemillas, r−r·cos(π s)(3.7) Para obtener el error medio por cada voto procesado, o lo que es lo mismo, la altura media del área subrayada en la Figura 3.10, el área del error tiene que ser dividida por el tamaño del lado 2y: Error medio por voto =área de error 2·y= π·r2 s−(r2·sin(π s)·cos(π s)) 2·r·sin(π s) =r· π s−(sin(π s)·cos(π s)) 2·sin(π s)=r·(π 2·s·sin(π s)−cos(π s) 2) =r 2·(π s·sin(π s)−cos(π s)) (3.8) 46 Capítulo 3. Detección de Patrones 3.6.2. Cota del error Como muestra la Figura 3.9, el número de semillas aplicadas en el rasterizador juega un papel fundamental para reducir el error en la interpolación de votos a través de líneas rectas (polígonos circunscritos dentro de la circunferencia). Este comportamiento se define mejor en la gráfica que acompaña a la Tabla 3.4, donde se representa el error máximo por voto frente al número de semillas, mostrando un decremento exponencial conforme aumenta el radio del círculo. Un umbral interesante aparece cuando el error máximo supera el valor correspondiente a la distancia de un píxel (ver Tabla 3.4 - línea horizontal sobre la gráfica). Por encima de este punto, no merece la pena incrementar el número de semillas porque el proceso de interpolación está enfocado hacia una representación en pantalla con la limitación del tamaño de píxel real. Puesto que un aumento del número de semillas repercute en la carga computacional de la GPU sin lograr unos resultados más precisos, este punto representa el valor óptimo para scuando se aplica el rasterizador. Dicho valor, acotado por un error máximo de un píxel, se obtiene mediante la siguiente expresión: r−r·cos(π s)<1⇒(1 −cos(π s)) ·r < 1⇒cos(π s)>r−1 s(3.9) La Tabla 3.4 muestra el número de semillas que cumplen la condición anterior para diferentes tamaños de círculo. El número de semillas empleadas por el rasterizador es siempre lo suficientemente bajo para garantizar un gran rendimiento computacional de la transformada Hough en términos de tiempo de ejecución, tal y como se verá en la Sección 3.9. 3.7. Optimizaciones en GPU Una desventaja de la implementación básica usando texturas bidimensionales para almacenar los votos es que todos los votos coincidentes que proceden de diferentes círculos o polígonos cuentan como un único voto. Este suceso es una consecuencia del empleo de primitivas básicas dentro de la misma pasada de renderizado. Para evitar esta pérdida de votos, el proceso de renderizado debe descomponerse en varias fases de renderizado, donde cada una de ellas sólo contiene puntos de contornos que se van a solapar en el espacio de parámetros. De esta forma, los spolígonos que representan el espacio de votado se renderizan alineados en una malla 2D, estando sus centros separados por 2r+1 píxeles, tanto en la componente horizontal como vertical (ver Figura 3.11a). Esta implementación hace uso de la primitiva GL_LINE_LOOP y 3.7. Optimizaciones en GPU 47 Tamaño radio Min. semillas r = 20 píxeles s = 10 r = 50 píxeles s = 16 r = 100 píxeles s = 23 r = 200 píxeles s = 32 Tabla 3.4: Número mínimo de semillas, s, para mantener el error máximo por debajo de la distancia de un píxel para diferentes tamaños de radio. La gráfica de la derecha representa el error en píxeles cuando se incrementa spara todos los radios. precisa de (2∗r+1)2pasadas de renderizado, consiguiendo un rendimiento mediocre para radios de gran tamaño. Estrategia Número de pasadas Estimación Precisión de la GPU implementada de renderizado del rendimiento GeForce 7 GeForce 8 Caso base (2r+ 1)2Deficiente Excelente Excelente Agrupando semicírculos 2(2r+ 1) Razonable Excelente Excelente Juntar canales de color 2(2r+1) 3Bueno Bueno Excelente Mezclado de votos 1 Excelente Deficiente Razonable Tabla 3.5: Conjunto de optimizaciones empleadas con el rasterizador de la GPU en términos de precisión y eficiencia. La Tabla 3.5 muestra la relación entre los votos perdidos y el factor de aceleración. El número de pasadas de renderizado de la segunda columna representa el peor caso, ya que en alguna de las fases de renderizado es posible que no exista solapamiento entre los diferentes puntos de contorno. Esta reducción será mayor para imágenes con pocos puntos, y podrían aplicarse una serie de optimizaciones en GPU con la aplicación previa de una fase de preprocesamiento sobre los contornos. Sin embargo, este capítulo estudia las vías de optimización para el peor caso como fiel indicador de la complejidad computacional de la implementación GPU. 48 Capítulo 3. Detección de Patrones (a) GL_LINE_LOOP (s = 8). (b) Agrupando semicírculos (primera mitad). (c) Agrupando semicírculos (segunda mitad). Figura 3.11: Planos de renderizado empleados en las distintas optimizaciones en GPU. (a) Polígonos de ocho caras alineados en una malla 2D cuyos centros están separados en 2r+ 1 píxeles para evitar la colisión de votos. (b) y (c) Todos los círculos “C” se renderizan en una primera pasada, mientras que los círculos “D” lo hacen en una segunda pasada adicional. 3.7. Optimizaciones en GPU 49 3.7.1. Agrupando semicírculos La gran cantidad de pasadas de renderizado que requiere la implementación básica puede ser reducida a 2·(2r+ 1) si los círculos son divididos en semicírculos y son procesados como se muestra en las Figuras 3.11b y 3.11c sobre el eje horizontal. Todos los semicírculos “C"son procesados en una primera pasada de renderizado y, posteriormente, se procede de igual manera con los semicírculos “D". Aún así, el eje vertical sigue requiriendo 2r+ 1 pasadas, pero ahora cada plano de renderizado acumula más votos y la implementación es más eficiente en GPU. La primitiva de interpolación empleada en esta implementación y en las sucesivas que de ella se derivan es GL_LINES_STRIP. 3.7.2. Uso de la caché de vértices y la rotación de matrices Los procesadores de vértices pueden beneficiarse de los datos procesados en las pasadas previas de renderizado si hacen uso de la caché de vértices y de la rotación de matrices para explotar las simetrías y la reutilización de los datos. Un ejemplo se muestra en la Figura 3.11c, que puede obtenerse a través de una rotación de todos los votos en un ángulo de 180 grados respecto a la pasada de renderizado anterior mostrada en la Figura 3.11b. De este modo, se evita toda la generación de vértices en GPU al igual que su comunicación con la CPU. La implementación de esta idea se lleva a cabo creando, en primer lugar, una lista de vértices en la CPU y, posteriormente, renderizando los vértices con DrawElements en lugar de mediante los bloques típicos GL_BEGIN-GL_END. La reorganización de los vértices consume tal cantidad de tiempo que no es posible amortizarla con la contribución de la caché de vértices debido a su reducido tamaño. 3.7.3. Mezcla de los canales de color Las pasadas de renderizado pueden reducirse a la tercera parte si se aprovechan los diferentes canales de color RGB para proyectar diferentes planos y finalmente fusionarlos en uno sólo. En este caso, la colisión de votos se evita, los procesadores de píxeles realizan la acumulación de votos para cada canal y la CPU se encarga de realizar la fusión de los tres canales. El resultado es un renderizado triple con una ejecución más rápida en GPU. 56 Capítulo 3. Detección de Patrones s. El decremento de δt implica un incremento proporcional en el número de votos a procesar por punto de contorno, y en este escenario la GPU muestra un procesamiento mucho más rápido, especialmente en imágenes que contienen pocos círculos y radios pequeños. Por otro lado, el decremento de stiene dos consecuencias para los círculos que son aproximados por polígonos con pocas caras: se interpolan más votos, sobrecargando el rasterizador, y se procesan menos votos, reduciendo la carga del procesador de vértices. En cargas de trabajo pequeñas (imágenes con 20 círculos), el decremento de spasa desapercibido en el tiempo de ejecución, lo que sugiere un rasterizado infrautilizado. En cargas de trabajo intensas (subconjunto de imágenes que contienen 500 círculos), el tiempo de ejecución se alivia alrededor de 1/3 cuando el número de semillas es la mitad. Votos emitidos por GL_POINTS Radio Imágenes representativas del conjunto de datos de entrada δt = 0.01 δt = 0.02 δt = 0.05 δt = 0.1 20 c20r20, c100r20, c500r20 628 314 125 63 50 c20r50, c100r50, c500r50 628 314 125 63 100 c20r100, c100r100, c500r100 628 314 125 63 Votos emitidos e interpolados por el rasterizador Radio Imágenes representativas del s = 4 s = 8 s = 16 s = 32 conjunto de datos de entrada em. in. em. in. em. in. em. in. 20 c20r20, c100r20, c500r20 4 109 8 114 16 108 32 93 50 c20r50, c100r50, c500r50 4 287 8 298 16 296 32 281 100 c20r100, c100r100, c500r100 4 561 8 604 16 608 32 595 Tabla 3.9: Número de votos emitidos (em.) e interpolados (in.) dependiendo de la implementación en GPU y el radio de los círculos (en píxeles) que va a ser detectado por la transformada de Hough. 3.10. Conclusiones El cauce de segmentación gráfico de una GPU puede aprovecharse en toda su extensión según la naturaleza del problema a resolver. Este capítulo presenta una implementación alternativa de la transformada de Hough aprovechando el conjunto de la tarjeta gráfica como una extraordinaria combinación de bajo coste y gran rendimiento. Nuestras implementaciones y optimizaciones de la THC contribuyen con un recorrido sobre las diferentes unidades funcionales de la GPU, desde una versión básica que sólo aprovecha el procesador de vértices para procesar los votos y el procesador de píxeles para acumularlos, hasta versiones más avanzadas donde (1) el rasterizador entra en juego para interpolar el espectro completo de votos, y (2) las unidades 3.10. Conclusiones 57 Casos de optimización Hardware Imagen CPU GPU, GL_POINTS GL_LINE_LOOP Tiempo Tiempo Aceleración Tiempo Aceleración c20r20 528.7 64.7 8.17x 38.8 13.62x c100r20 745.8 89.0 8.37x 50.8 14.68x c500r20 1770.6 23.5 7.51x 117.3 15.09x c20r50 1535.3 130.0 11.81x 44.4 34.57x c100r50 2871.0 303.0 9.47x 78.2 36.71x c500r50 9002.1 1084.0 8.30x 237.4 37.91 c20r100 3620.9 293.0 12.35x 53.5 67.68x c100r100 6822.4 684.0 9.97x 94.8 71.96x c500r100 29190.5 3737.0 7.81x 396.8 73.56x Versiones con técnicas de optimización Imagen Agrupando semicírculos Mezcla de canales Blending de votos Tiempo Aceleración Tiempo Aceleración Tiempo Aceleración c20r20 19.2 27.53x 11.8 44.80x 2.0 264.35x c100r20 24.6 30.31x 15.1 49.39x 9.8 76.10x c500r20 69.4 25.51x 58.6 30.21x 48.2 36.73x c20r50 46.0 33.37x 28.5 53.87x 5.0 307.06x c100r50 62.9 45.64x 37.2 77.17x 24.8 115.76x c500r50 164.4 54.75x 143.6 62.68x 117.7 76.48x c20r100 92.3 39.22x 57.6 62.86x 10.1 358.50x c100r100 108.2 62.82x 67.1 101.67x 34.4 198.32x c500r100 304.2 95.95x 254.9 114.51x 213.8 136.53x Tabla 3.10: Tiempos de ejecución (en milisegundos) y factores de aceleración de todas las optimizaciones desarrolladas para la transformada de Hough. Las versiones de CPU y GL_POINTS se han ejecutado para un ángulo δt tal que la cantidad de votos emitidos es similar a las versiones con rasterizador. Los tiempos del rasterizador se muestran para un valor s de 32 semillas. Los factores de aceleración mínimos se han subrayado y los máximos se han resaltado. de blending se utilizan para almacenar los votos de cada textura de una manera más eficiente. Para la versión que emplea el rasterizador del procesador de vértices haciendo uso de la primitiva GL_LINES, hemos desarrollado una fórmula para determinar el mínimo número de semillas que garantiza que el error máximo siempre estará por debajo del umbral de la distancia de un píxel y que nos permite ahorrar en tiempo de ejecución sin deterioro en la calidad de los resultados. La fórmula es una función de dos parámetros que influyen negativamente en la precisión del proceso de votación: tamaño superior del radio para la detección de círculos, y resolución inferior para las texturas usadas en la acumulación de votos en GPU. En una cantidad similar de votos, el conjunto completo de optimizaciones en GPU 58 Capítulo 3. Detección de Patrones Imagen Tiempo CPU GL_POINTS GL_LINE_LOOP Blending c20r20 266.7 44.6 29.9 (4) 1.9 (4) (δt = 0.1) 30.7 (8) 1.9 (8) 528.7 64.7 34.1 (16) 2.0 (16) (δt = 0.05) 38.8 (32) 2.0 (32) c100r20 375.3 63.2 32.4 (4) 8.0 (4) (δt = 0.1) 35.9 (8) 8.0 (8) 745.8 89.0 39.9 (16) 8.5 (16) (δt = 0.05) 50.8 (32) 9.8 (32) c500r20 888.2 14.0 45.4 (4) 17.2 (4) (δt = 0.1) 54.9 (8) 23.2 (8) 1770.6 23.5 74.0 (16) 32.7 (16) (δt = 0.05) 117.3 (32) 48.2 (32) c20r50 303.9 51.2 31.6 (4) 4.8 (4) (δt = 0.1) 33.4 (8) 4.8 (8) 1535.3 130.0 35.5 (16) 4.9 (16) (δt = 0.02) 44.4 (32) 5.0 (32) c100r50 565.7 80.0 37.0 (4) 20.2 (4) (δt = 0.1) 46.4 (8) 20.3 (8) 2871.0 303.0 56.3 (16) 21.6 (16) (δt = 0.02) 78.2 (32) 24.8 (32) c500r50 1771.3 256.0 64.7 (4) 43.1 (4) (δt = 0.1) 81.4 (8) 54.7 (8) 9002.1 1084.0 133.9 (16) 77.9 (16) (δt = 0.02) 237.4 (32) 117.7 (32) c20r100 359.8 59.9 34.5 (4) 9.7 (4) (δt = 0.1) 35.1 (8) 9.7 (8) 3620.9 293.0 40.0 (16) 9.9 (16) (δt = 0.01) 53.5 (32) 10.1 (32) c100r100 680.6 97.0 39.2 (4) 30.9 (4) (δt = 0.1) 47.7 (8) 31.0 (8) 6822.4 684.0 61.5 (16) 32.0 (16) (δt = 0.01) 94.8 (32) 34.4 (32) c500r20 2919.0 403.0 87.8 (4) 105.6 (4) (δt = 0.1) 121.3 (8) 109.6 (8) 29190.5 3737.0 216.2 (16) 145.7 (16) (δt = 0.01) 396.8 (32) 213.8 (32) Tabla 3.11: Tiempos de ejecución (en milisegundos) para diferentes estrategias y carga de trabajo de la THC: los tiempos de CPU y GL_POINTS evalúan dos pasos de ángulo diferentes, δt, y GL_LINE_LOOP y Blending evalúan diferentes cantidades de semillas (proporcionadas entre paréntesis). Los valores con un número similar de votos procesados son resaltados para cada imagen con el objetivo de hacer una comparación lo más justa posible entre estrategias. 3.10. Conclusiones 59 desarrollado en este capítulo acelera el rendimiento de CPU en un factor entre 25x y 358x, dependiendo del tamaño del radio y el número de círculos que contienen las imágenes de muestra. Finalmente, con el objetivo de realizar una comparativa entre CPU y GPU en términos de rendimiento y escalabilidad, hemos llevado a cabo una serie de optimizaciones en CPU. Con la ayuda de tablas trigonométricas precalculadas, hemos obtenido un factor de mejora de 4x para ambas versiones. Si empleamos técnicas de gradiente, el rendimiento obtiene una ganancia de 40x a costa de aumentar la sensibilidad al ruido. En este aspecto, la GPU cosecha ejecuciones más rápidas, y siempre con la robustez del algoritmo clásico de la THC. La escalabilidad se ha analizado sobre la misma familia de plataformas hardware entre 2006 y 2008. La GPU GeForce 8800 obtiene un rendimiento entre 10 y 20 veces superior a la GPU GeForce 7950 GX2, mientras que la CPU Intel Core 2 Duo reduce los tiempos de ejecución en aproximadamente la mitad respecto al Pentium 4. Parte III GPGPU contemporánea 61 4Clasificación celular de los tumores neuroblásticos 4.1. Análisis de imagen histopatológico El neuroblastoma es un cáncer agresivo del sistema nervioso que afecta principalmente a los niños. El pronóstico de la enfermedad depende en gran medida del examen histopatológico de los tejidos seccionados y tintados con hematoxilina y eoxina (H & E), realizado visualmente con un microscopio por patólogos expertos. El sistema de clasificación actual se basa en algunas características morfológicas de la histología del tejido tales como el grado de desarrollo del stroma Schwannian, el grado de diferenciación y la mitosis y el índice karryorhexis. Dicho análisis junto con otros indicadores como los factores genéticos y la edad del paciente puede desembocar en un pronóstico favorable o desfavorable. Sin embargo, el examen cuantitativo de las muestras de tejido depende frecuentemente de variaciones entre los expertos y de un sesgo muestral. Por tanto, dicho método propenso a errores puede ocasionar la planificación de un tratamiento inadecuado para los pacientes. Para solucionar este inconveniente, existen métodos de análisis cuantitativos por ordenador que son más objetivos y precisos a la hora de extraer las características diagnósticas que desde el punto de vista del ser humano sería imposible apreciar [13, 27, 70]. Este capítulo va a tratar la clasificación del grado de desarrollo del stroma de Schwannian, uno de los indicadores más importantes de malignidad según el Sistema Internacional de Clasificación del Neuroblastoma (INCS). El INCS requiere distinguir entre estroma pobre y rico de los tejidos, lo que vamos a abordar como un problema 63 64 Capítulo 4. Clasificación celular de los tumores neuroblásticos de clasificación de patrones de los tejidos. La Figura 4.1 muestra regiones pequeñas de tejido rico y pobre en estroma. La textura del tabique del estroma (ver Figura 4.1a) está bastante diferenciada respecto a los neutrófilos (ver Figura 4.1b). Las estructuras de fibrina del estroma con forma de pelo muestran patrones organizados localmente en direcciones concretas, mientras que el neutrófilo sigue una estructura de malla. Esta información se ha extraído de las características de la textura a través de estadísticas de co-ocurrencia y patrones binarios locales (Local Binary Patter, LBP). Una descripción más detallada de todo el análisis de la imagen para la clasificación del desarrollo del estroma puede encontrarse en [70]. (a) Tejido rico en estroma. (b) Tejido pobre en estroma. Figura 4.1: Muestras de imágenes de neuroblastoma. Debido a las grandes resoluciones de las imágenes procesadas (hasta más de 100K x 100K, con tamaños superiores a 30 gigabytes), el conjunto completo de imágenes seccionadas se divide en grupos más pequeños que no se solapan, y al que se le aplica el análisis de imagen rutinario de forma independiente. Posteriormente, los resultados de cada imagen se fusionan para obtener el mapa de clasificación correspondiente, que consiste en un conjunto de etiquetas de clasificación para cada imagen (estroma pobre o rico). La resolución del mapa de clasificación dependerá de la resolución seleccionada para los grupos de imágenes resultantes de la visión, de forma que las imágenes de menor resolución originan un mapa de clasificación de mayor resolución, y viceversa. La Sección 4.1.1 resume el análisis de imagen propuesto para cada imagen. 4.1.1. El algoritmo Tal y como ya se ha mencionado, el objetivo es la diferenciación en las texturas entre el tejido rico en estroma y el tejido pobre en estroma. Para ello se ha creado un vector de características compuesto por estadísticos de co-ocurrencia y LBPs para 4.1. Análisis de imagen histopatológico 65 cada imagen. La clasificación en este espacio de características se finaliza de forma supervisada. La Figura 4.2 presenta el flujo de trabajo del algoritmo de análisis de imagen para la clasificación del desarrollo del estroma en imágenes de neuroblastoma tintadas con hematoxilina y eoxina. Dicho análisis se compone de cuatro fases que se definen a continuación. Figura 4.2: Diagrama de flujo del algoritmo de clasificación del estroma. Transformación del espacio de colores. Para obtener una mejor representación de la información de la intensidad y de los colores de manera independiente, el espacio de colores LA*B* proporciona un espacio de colores perceptiblemente uniforme. La uniformidad perceptible implica que un cambio de color de la misma proporción siempre se manifiesta visualmente con la misma importancia [60]. En el espacio de colores LA*B*, el canal L corresponde a la luminancia o iluminación y los canales A* y B* corresponden a dos dimensiones opuestas de color (los colores opuestos son aquellos que son percibidos de esta manera por el ojo humano). La información de iluminación y de colores se separa para que el procesamiento de la textura sea más apropiado para el reconocimiento de patrones a posteriori. Características de co-ocurrencia. Los estadísticos de co-ocurrencia obtenidos de las matrices de co-ocurrencia facilitan cuatro características por canal: contraste, correlación, homogeneidad y energía. Dichas características se tipifican como las características de Haralick [77]. Las matrices de co-ocurrencia miden la frecuencia con la que un píxel con una intensidad iaparece en una relación espacial específica con un píxel de intensidad j(ver Figura 4.3). 72 Capítulo 4. Clasificación celular de los tumores neuroblásticos un papel fundamental para la eficiencia de la GPU, mientras que la influencia en la CPU es marginal. Aunque la CPU obtiene resultados más rápidos con matrices de 16 x 16, el empleo de matrices de 8 x 8 ó 4 x 4 conducen a unos resultados de clasificación similares (ver Tabla 4.1) y la computación en GPU es mucho más rápida. En general, el tiempo de ejecución para la aplicación del análisis de las imágenes patológicas se puede reducir en un factor 321x (matriz de co-ocurrencia 4 x 4 en Tabla 4.3), donde un factor de 7x proviene de la traslación del código Matlab a C++ dentro de la CPU y el resto corresponde a la contribución de la GPU. Mientras que la primera mejora se mantiene constante, la segunda es sensible al tamaño de la matriz de co-ocurrencia y la forma en que se procesa en la GPU. Para simplificar el estudio, de aquí en adelante todos los resultados se refieren a una matriz de co-ocurrencia de 4 x 4, ya que se obtiene la mejor aceleración sin perder precisión en los resultados de clasificación. Tamaño de la matriz de co-ocurrencia Matlab/C++ C++/GPU Matlab/GPU 16 x 16 7.1x 5.1x 36.4x 8 x 8 7.1x 17.2x 122.6x 4 x 4 7.0x 45.3x 321.2x Tabla 4.3: Factor de aceleración para diferentes matrices de co-ocurrencia sobre una imagen de tamaño 1020 x 916 píxeles bajo diferentes plataformas. Las Tablas 4.4 y 4.5 y la Figura 4.5 muestran cómo la mejora en CPU no se mantiene constante para cada tarea cuando el tamaño de la imagen va incrementándose: las ganancias se sustentan en el cálculo de las características estadísticas, donde LA*B* y LBP reducen la ganancia como penalización por interpretar el código Matlab en tiempo de ejecución, hecho éste que se reduce cuanto más datos por instrucción surgen del programa. Además, el código Matlab se comporta mejor trabajando con la memoria de forma masiva por el buen uso de la memoria virtual, mientras que el código C++ resulta más efectivo en imágenes pequeñas, haciendo un uso más eficiente de la memoria caché. Las Figuras 4.6 y 4.7 muestran la ganancia adicional de la GPU sobre el código Matlab y C++. El cómputo de la conversión de color LA*B* es de nuevo la tarea con un impacto de rendimiento mayor, ya que constituye un operador streaming típico mientras retiene la mayor parte de la carga de trabajo (ver Tabla 4.5). Este papel de liderazgo se mantiene también en la ganancia total, aunque las características estadísticas no rindan mejor en GPU que en C++ para tamaños de imagen muy pequeños, donde el cómputo no amortiza el coste de mover los datos hacia la CPU. Por último, los factores de aceleración se mantienen para imágenes grandes (1 x 1 K y 2 x 2 K). Finalmente, la Tabla 4.4 muestra hasta qué punto el peso computacional de cada 4.2. Implementación pre-CUDA en GPU: Programación gráfica clásica 73 Etapa de En la CPU: Entre procesadores: Global: procesamiento C++/Matlab CPU/GPU Matlab/GPU Conversión LA*B* 5.9x - 5.2x 69.2x - 1409.6x 406.1x - 7391.3x Estadísticos 122.2x - 90.0x 0.2x - 2.1x 21.8x - 192.1x LBP 8.3x - 3.9x 4.2x - 38.3x 34.6x - 148.3x Total 13.3x - 7.6x 2.6x - 46.3x 33.4x - 350.9x Tabla 4.4: Influencia de la carga de trabajo en las ganancias de rendimiento para cada una de las tareas involucradas en el análisis del neuroblastoma. Las ganancias se muestran para imágenes de 128 x 128. El tamaño de las matrices de co-ocurrencia es de 4 x 4 en todos los casos. Equipamiento hardware del año 2006 Desaceleración relativa ( %) CPU Mobile Intel Core Duo @ 2.16GHz 45.2 CPU Desktop Intel Pentium 4 @ 3.4GHz 34.1 CPU Desktop Mac Pro Intel Xeon @ 2.66GHz 12.0 GPU Nvidia GeForce 96.5 Tabla 4.5: Comparativa de la escalabilidad del hardware desde el lado de la CPU y la GPU. La última columna resume la desaceleración en comparación con el hardware del año 2007 mostrado en la Tabla 4.2. El tamaño de la imagen es de 1020 x 960 píxeles. Figura 4.5: Factores de aceleración entre la versión C++ y Matlab para diferentes tamaños de imagen (en píxeles). tarea depende de la implementación real. En la CPU, la carga principal viene dada por la conversión del espacio de color, mientras que en la GPU las características estadísticas componen la mayoría del tiempo de ejecución. Como se muestra en la última columna, la implementación en GPU proporciona una aceleración de más de dos órdenes de magnitud para imágenes de 2 x 2 K. 74 Capítulo 4. Clasificación celular de los tumores neuroblásticos Figura 4.6: Factores de aceleración entre la versión GPU y C++ para diferentes tamaños de imagen (en píxeles). Figura 4.7: Factores de aceleración entre la versión GPU y Matlab para diferentes tamaños de imagen (en píxeles). Variación del rendimiento en función del hardware. En CPU se ha ejecutado el código en tres equipos diferentes del año 2006: una CPU portátil, un Pentium 4 convencional y un Xeon de gama alta. Por parte de la GPU, el código se ha ejecutado en una GeForce 7950 GX2 en representación del hardware gráfico del 2006. Los resultados se muestran en la Tabla 4.5, donde se puede apreciar que la mejora de rendimiento de la GPU es mucho mayor que en todas las CPUs. Mientras que el modelo de GPU del año 2007 consigue un factor de ganancia sobre 2x en comparación con el modelo de 2006, la CPU del 2007 obtiene apenas una aceleración de un 1.1x frente a un modelo de gama alta del 2006, y sólo un 1.3x aproximadamente sobre la versión portátil del 2006. Validación numérica. Para juzgar la similitud de los valores procesados en diferentes plataformas, ejecutamos el código de análisis de imágenes en cada plataforma, 4.2. Implementación pre-CUDA en GPU: Programación gráfica clásica 75 obteniendo un vector de características para 500 imágenes de nuestro conjunto de imágenes de entrenamiento. La Tabla 4.6 muestra el valor medio y la desviación estándar de las diferencias entre los vectores de 13 características sobre las 500 imágenes. Las diferencias están alrededor de un 3 % con un valor menor de 0.02 en promedio, siendo estas pequeñas variaciones inapreciables en la precisión de la clasificación final (ver Tabla 4.1). Versiones comparadas Diferencia media Desviación estándar Media del error relativo máximo Matlab/C++ 0.00014 - 0.012 0.00018 - 0.01 2.00 C++/GPU 0.00065 - 0.021 0.00043 - 0.05 3.46 Matlab/GPU 0.0015 - 0.017 0.00075 - 0.05 Tabla 4.6: Precisión de los valores de salida para diferentes plataformas hardware. El rango de variación corresponde a un test sobre 500 imágenes para las 12 características estadísticas junto con el operador LBP (13 valores en total). Manejo de la memoria y tamaño de las imágenes. La capacidad de la memoria DRAM se ha duplicado cada año y medio, en línea con la predicción de la Ley de Moore. Con el tamaño actual de 1 GB de memoria de vídeo y el formato de color SRGB de 8 bits usado en la GPU para representar las imágenes como texturas, nuestra aplicación se podría ejecutar completamente en un solo núcleo del año 2012. Hasta entonces había que emplear alguna lógica de particionado porque la GPU no utilizaba memoria virtual. Un particionado en imágenes más pequeñas puede beneficiarse de las cachés de la GPU de manera similar a cómo los compiladores aplican las transformaciones de bloques en el ámbito de la CPU, con la distinción que la memoria caché de la GPU es mucho más pequeña y la memoria DRAM mucho más rápida. La Tabla 4.7 muestra los factores de aceleración conseguidos usando diferentes tamaños para el particionado y para una matriz de co-ocurrencia de tamaño 4 x 4 para la extracción de las características estadísticas. El ancho de banda determina un particionado máximo en imágenes de 128 x 128 píxeles, lo que supone el límite crítico donde el tiempo de inicialización se amortiza por la velocidad de la comunicación. Por otro lado, la capacidad de memoria y la programación del software determinan el tamaño máximo de las imágenes, que en la práctica fue de 8 x 8 K para la GPU usando OpenGL 2.1.4, 4 x 4 K para la CPU usando C++, y 2 x 2 K cuando se usa Matlab de 32 bits. La Figura 4.7 muestra que el rendimiento mejora con el tamaño de las imágenes, reflejando que los tamaños de 1 x 1 K y 2 x 2 K son prioritarios de cara al factor de aceleración. Además, se ha estudiado el efecto de un particionado en imágenes no cuadradas. Cuando la dimensión horizontal de las imágenes es mayor, la CPU es 2.7 % más rápida usando Matlab, lo que refleja que el índice de los arrays sigue un orden basado en las 76 Capítulo 4. Clasificación celular de los tumores neuroblásticos Tamaño de imagen Matlab/C++ C++/GPU Matlab/GPU 128 x 128 13.3x 2.6x 33.4x 256 x 256 7.0x 9.1x 64.2x 512 x 512 7.5x 24.1x 182.1x 1024 x 1024 7.6x 439.4x 297.7x 2048 x 2048 7.6x 46.3x 350.9x Tabla 4.7: Factores de aceleración obtenidos para una matriz de co-ocurrencia de 4x4 sobre diferentes tamaños de imagen (expresados en píxeles). columnas en la implementación de Matlab que se beneficia de un patrón de localidad a la hora de realizar las peticiones de las líneas de caché. Sin embargo, la GPU no obtiene ningún beneficio usando imágenes no cuadradas. Comparativa global. Para componer el cómputo de todas las imágenes de la muestra inicial, el tiempo total de ejecución para la clasificación de la región del estroma para una imagen de 50 x 50 K alcanza las 4 horas y 45 minutos en Matlab, 37 minutos para la ejecución C++ en CPU y sólo 145 segundos para la implementación más rápida en GPU. Considerando que aproximadamente 600 pacientes son diagnosticados con neuroblastoma al año y que cada tumor lleva consigo aproximadamente 5-6 biopsias recogidas del paciente, el tiempo total de procesamiento de estas imágenes sería de 21 meses usando Matlab en una CPU de doble núcleo frente a unos 5.3 días para la versión GPU descrita en este capítulo. La Sección 4.2.1 resume la mejor estrategia de distribución de la carga entre los procesadores: la GPU procesa la conversión LA*B* y el operador LBP, mientras que la CPU procesa las sumas de las 12 características estadísticas. En una CPU de doble núcleo programada con POSIX threads, uno de los núcleos puede trabajar en el algoritmo y en interactuar con la GPU, mientras el otro núcleo se encarga de cargar la siguiente imagen a procesar. De esta forma se ocultaría la latencia de entrada/salida para el proceso completo. En general, una vez que la conversión LA*B* ha terminado, el operador LBP y las características estadísticas son independientes y pueden ser procesadas en paralelo en la GPU y el primer núcleo de la CPU, respectivamente. Mientras tanto, el segundo núcleo puede encargarse de cargar el resto de las imágenes. 4.3. Implementación del código de análisis de imagen en CUDA La implementación CUDA se ha llevado a cabo con su ciclo de desarrollo típico. En primer lugar, el código se compila con las opciones de compilación que generan un informe sobre el aprovechamiento de los recursos hardware de cada kernel (regis- 4.3. Implementación del código de análisis de imagen en CUDA 77 tros requeridos por cada hilo y memoria compartida ocupada por cada bloque). Dicho informe analizado en detalle proporciona el número de hilos y bloques que son necesarios lanzar en la ejecución para conseguir la máxima ocupación dentro de cada multiprocesador. Sí aún así no hay posibilidad de obtener una eficiencia satisfactoria, el código debe revisarse para reducir la presión sobre la memoria. Debido al gran rendimiento que la GPU proporciona con operaciones en punto flotante, los accesos a memoria suelen ser los cuellos de botella del rendimiento en diversas partes de la aplicación. La imagen de entrada (1K x 1K x 3 bytes) es bastante más grande que el tamaño de la memoria compartida (16 KB), así que la prioridad está en las estructuras de datos como las matrices de co-ocurrencia (fase 2) y los histogramas parciales (fase 4). Sin embargo, en la fase 1, aunque el procesamiento es un barrido sobre todos los píxeles, el tiempo de ejecución es menor cuando se usa la memoria compartida (2.32 ms contra 2.77 ms - ver Tabla 4.8). Además, en la fase 3, los píxeles de entrada se trasladan a memoria compartida porque el cálculo del operador LBP reutiliza mucho los datos. Para ilustrar la progresión de una implementación CUDA de una manera didáctica y detallada, se van a ilustrar las diferentes optimizaciones que hemos llevado a cabo en cada fase. 4.3.1. Detalles de implementación en CUDA Fase 1: Conversión de color Como punto de partida se emplean tipos de datos float3 de 24 bits para cada canal de color. Sin embargo, usando un mecanismo de relleno para que el ancho de datos sea de 32 bits se puede mejorar el rendimiento en todas las optimizaciones en las que entra en escena la memoria compartida (ya que la salida de un banco de esta memoria es de 32 bits). La intención es eliminar la famosa coalescencia de datos, que para esta fase permite ahorrar un 35% del tiempo de computación, aunque a costa de aumentar el tiempo de comunicación (ver Tabla 4.8). A continuación puede observarse cómo un tipo de datos uchar de 8 bits es suficiente para la precisión de la aplicación. Tal y como se espera, el tiempo de comunicación se reduce aproximadamente a la cuarta parte. El uso de las opciones de compilación que informan sobre el uso de los registros y la memoria en CUDA indican que el kernel CUDA de la conversión de color hace uso de 13 registros y 1064 bytes de memoria compartida, permitiendo una ocupación máxima del 75 % del multiprocesador para una ejecución en la que el bloque tenga entre 176 y 192 hilos. Sin embargo, tomamos la decisión de seleccionar 256 hilos, lo 78 Capítulo 4. Clasificación celular de los tumores neuroblásticos Fase Descripción / Optimizaciones Tiempo de ejecución (ms.) Comun. Comput. Total 1: Versión base: un float3 para cada canal de color 8.49 3.71 12.20 Conversión Coalescencia (insertando el canal alfa) con float3 10.79 2.44 13.23 de color Remplazo de float3 por uchar (256 hilos/bloque) 2.98 2.77 5.75 RGB a Uso de memoria compartida (256 hilos/bloque) 2.98 2.32 5.30 LA*B* Utilizando entre 169 y 192 hilos/bloque 2.98 2.43 5.41 2: Versión base: memoria global 15.40 15.40 Parámetros Memoria compartida para las matrices de co-ocurrencia 4.48 4.48 estadísticos Conflictos resueltos en bancos de mem. compartida 2.58 2.58 3: Versión base: Hilos especiales sobre los bordes (halos) 2.29 2.29 Operador Bloques de 16 x 16 hilos, computación en 14 x 14 1.82 1.82 LBP Bloques de 8 x 8 hilos, computación en 6 x 6 2.31 2.31 4: Versión base: memoria global 4.02 2.08 6.10 Histograma Memoria compartida para histogramas locales 0.31 0.61 0.92 Conflictos entre warps resueltos en bancos de mem. 0.31 0.59 0.90 Tiempo Versión base 12.51 13.48 35.99 total Con memoria compartida 3.29 9.23 12.52 GPU Versión óptima 3.29 7.31 10.60 Tiempo total CPU: (1:) 880.31 ms + (2:) 43.24 ms + (3:) 156.28 ms + (4:) 7.49 ms = 1087.32 ms. Tabla 4.8: Principales optimizaciones CUDA en la aplicación de análisis de imagen. Los tiempos de ejecución corresponden a una solo imagen de 1K x 1K píxeles en una única GPU. La fase 1 muestra los tiempos para la comunicación entre CPU y CPU, y la computación real en GPU. Las fases 2 y 3 sólo computación (no existe comunicación). La fase 4 muestra los tiempos para la comunicación entre GPU y CPU después de la computación real en GPU. que supone un canje entre la ocupación (desde 75 % hasta 67 %) para obtener un mejor balanceo de la carga, ya que 256 es múltiplo de 32 (número máximo de hilos por warp), y a su vez es divisor de 768 (número máximo de hilos por multiprocesador) y 1024 (número máximo de píxeles por imagen). El resultado muestra que el tiempo de ejecución mejora ligeramente (ver Tabla 4.8). Desafortunadamente, el rendimiento máximo está limitado porque cada hilo necesita 11 registros, provocando que se usen menos de los 768 hilos de ocupación máxima por multiprocesador. El tiempo de ejecución óptimo para esta fase es de 2.98 ms. para la transferencia de los píxeles, y de 2.32 ms. para computar la conversión de color, tal y como se refleja en la Tabla 4.8. Fase 2: Características estadísticas Este kernel necesita 9 registros y 4132 bytes de memoria compartida, de manera que se pueden alojar hasta 3 bloques de 256 hilos en paralelo. Empleando 768 hilos por bloque podemos subir el factor de ocupación los recursos hardware de la GPU 4.3. Implementación del código de análisis de imagen en CUDA 79 hasta el 100%. Los píxeles se distribuyen equitativamente entre los hilos para procesar simultáneamente las matrices de co-ocurrencia. Finalmente, los resultados parciales se acumulan a través del operador de reducción. Dentro de esta fase, el resto consiste en procesar las matrices de co-ocurrencia evitando los conflictos de acceso a los 16 bancos de memoria. Con una matriz de 256 hilos distribuida en 16 x 16, el despliegue más sencillo de hilos sería el de forzar a los 32 hilos de un warp el acceso a sólo 8 bancos de memoria compartida, cuyo rendimiento estaría mermado por dicha limitación de paralelismo. Actuando de manera inteligente, el acceso de los hilos activos se puede distribuir de manera que ningún hilo tenga que esperar para acceder a su banco de memoria (ver Figura 4.8). Esta compleja optimización resuelve todos los conflictos en el acceso a memoria, reduciendo el tiempo de ejecución a 2.58 ms. desde 4.48 ms. Sin el uso de memoria compartida, una implementación base tardaría 15.40 ms (ver Tabla 4.8). Fase 3: Operador LBP La implementación del operador LBP conlleva la aplicación de una máscara de convolución de tamaño 3 x 3 píxeles seguida de una conversión de binario a decimal (ver Figura 4.4). Cada hilo necesita 10 registros y cada bloque de hilos usa 296 bytes de memoria compartida. Debido a las características del uso de la memoria, se puede alcanzar una ocupación de 100% con una distribución de 256 hilos desplegados en cuadrículas de 16 x 16. Cada hilo lee un píxel desde memoria global y lo almacena en la estructura de datos de memoria compartida. Los hilos situados en los bordes de las cuadrículas no pueden realizar procesamiento alguno puesto que no tienen acceso a todos sus vecinos más cercanos (ver Figura 4.9). El operador LBP para las regiones de los bordes se calcula con el siguiente bloque de hilos, para ello es necesario que haya solapamiento entre bloques de dos filas y dos columnas. Esta estrategia de cómputo sufre las consecuencias de un 23 % de ciclos ociosos y del correspondiente acceso a memoria redundante, a costa de conseguir un procesamiento más homogéneo y que los hilo estén menos sobrecargados. En comparación con la versión heterogénea, donde no existen ciclos ociosos y accesos redundantes a memoria, la versión homogénea es un 25% más rápida, llevando el mejor tiempo de ejecución hasta los 1.82 ms. Dado que el procesamiento del operador LBP es muy regular, existe la posibilidad de seleccionar diferentes tamaños para el despliegue de hilos, siempre y cuando la distribución se haga en cuadrículas simétricas. Por tanto, se ha investigado el efecto de las diferentes relaciones de hilos por bloque, en particular con los tamaños 20 x 20 (2.02 ms.), 18 x 18 (1.90 ms.), 16 x 16 (1.82 ms. - óptimo), 14 x 14 (1.89 ms.), 12 x 12 (2.00 ms.), 10 x 10 (2.16 ms.), y 8 x 8 (2.31 ms. - ver en Tabla 4.8). Un 80 Capítulo 4. Clasificación celular de los tumores neuroblásticos Figura 4.8: Matrices locales de co-ocurrencia en CUDA. La asignación de bancos en memoria compartida a cada hilo evita conflictos de memoria al procesar las matrices de co-ocurrencia. tamaño de cuadrícula de 16 x 16 es en la mayoría de los casos una buena elección para maximizar el rendimiento en CUDA. El estudio sobre los diferentes tamaños sólo pretende cuantificar la penalización que conlleva una mala elección de tamaño. Fase 4: Histograma La implementación de esta fase se basa en el núcleo del histograma incluido en la librería CUDA [64], donde la reducción global se delega a la CPU. Esta solución se ha considerado acertada porque el histograma sólo se computa una vez por imagen y no implica un gran tiempo de ejecución. El tiempo de ejecución para esta última fase es de 0.9 ms. La Tabla 4.8 resume todas las optimizaciones desarrolladas en cuanto a CUDA se 4.3. Implementación del código de análisis de imagen en CUDA 81 Figura 4.9: Operador LBP en CUDA. La memoria empleada por cada bloque es superior a la cantidad de hilos para poder computar las convoluciones en todos los bloques de forma homogénea. refiere, desglosando por fases y tiempos de ejecución. Dado que el modelo de programación CUDA para la versión empleada no permite la ejecución de kernels concurrentes entre multiprocesadores, las cuatro fases se ejecutan secuencialmente y la GPU se explota en su completitud para cada fase independientemente. En general, la aplicación de CUDA para esta aplicación de análisis de imágenes permite la reducción de los tiempos de ejecución en un factor entre 3 y 5 respecto a los tiempos obtenidos en la versión Cg (ver Tabla 4.9 para los tiempos sobre una imagen completa), con un factor adicional de 3x cuando se habilita la memoria compartida, añadiendo además un 20 % extra cuando se solventan los conflictos en los accesos a los bancos de memoria. 4.3.2. Datacutter La paralelización del análisis de imagen que hemos abordado en este capítulo ha contado con la ayuda del middleware Datacutter. Datacutter [6] ha sido desarrollado 88 Capítulo 4. Clasificación celular de los tumores neuroblásticos do desde la tecnología más clásica para programar las GPUs con la librería gráfica OpenGL y una plataforma híbrida sencilla de una única CPU y GPU, se empiezan a obtener claras evidencias de que el uso de la GPU para este tipo de análisis biomédicos permite ahorrar grandes cantidades de tiempo, aunque no por igual en todas las fases de computación. El diseño de la aplicación biomédica comienza en un lenguaje de alto nivel interpretado (Matlab) para definir un prototipo que posteriormente se traslada al lenguaje de programación C++ con uso exclusivo de CPU y otra versión con uso de GPU. Tras un análisis de las operaciones más adecuadas para cada tipo de procesador, el primer intento de plataforma híbrida consigue tiempos de ejecución alrededor de 45 veces más rápidos respecto a la versión C++ ejecutada exclusivamente en CPU, y 321x más rápidos frente a la versión prototipo en Matlab. La GPU es especialmente ventajosa cuando las operaciones a realizar requieren un gran ancho de banda, mientras que otras tareas como las matrices de co-ocurrencia son más apropiadas para la CPU. Un entorno más complejo compuesto de un clúster cooperativo de CPUs y GPU, y empleando la tecnología CUDA y DataCutter para paralelizar la computación internamente y entre los nodos, supone una sólida plataforma multiprocesador cooperativa y heterogénea en la que se puede explotar paralelismo a múltiples niveles. Concretamente, hasta cuatro niveles de granularidad de la arquitectura hardware y la aplicación software: (1) multi-nodo (usando DataCutter para el particionado de datos entre nodos), (2) SMP y a nivel de hilo (usando DataCutter para utilizar todos los recursos dentro de cada nodo y de cada procesador), (3) SIMD (usando CUDA para llenar todos los multiprocesadores de las GPU), y finalmente, (4) ILP (Paralelismo a nivel de instrucción), configurando los bloques de hilos dentro de la GPU de manera que siempre exista trabajo en la cola de procesos. Los resultados experimentales muestran el éxito de las técnicas empleadas, comenzando con el decremento de los tiempos de ejecución en un solo nodo CPU/GPU por el uso de diferentes optimizaciones intra-nodo. La ganancia se extiende al paralelismo inter-nodo para una ejecución escalable del sistema multiprocesador. En el análisis del conjunto de imágenes más grandes, para la configuración de clúster de 16 nodos, la implementación más simple de GPU DataCutter-CUDA es 31.3 veces más rápida que su implementación serie en CUDA. Al introducir dos GPUs por nodo, el tiempo de procesamiento de cada nodo cae por debajo del minuto, de manera que si se excluye el tiempo de entrada y salida y la descompresión de imágenes, queda demostrado su potencial. Además, cuando se emplea DataCutter para solapar los tiempos de computación con los de entrada y salida y descompresión, la GPU permanece con carga de trabajo constante. El resultado de esta aproximación llega a acelerar los tiempos en un factor de 12.94 sobre 16 nodos. 5Registro no rígido de imágenes microscópicas La caracterización de fenotipos asociados a genotipos específicos es fundamental para esclarecer las actividades de los genes y su interacción. La morfología de la estructura celular y tisular son aspectos del fenotipo que proporcionan la información necesaria para conocer procesos biológicos esenciales, tales como la iniciación del cáncer en el microentorno del tumor y la formación de redes neuronales en la corteza del cerebro. Sin embargo, las técnicas actuales que obtienen en la información tridimensional magnificada desde las muestras biológicas son bastante limitadas. Para este propósito suele emplearse un microscopio de fluorescencia multifotón y cofocal que necesita marcadores fluorescentes y tiene un campo de visión limitado. Por tanto, la reconstrucción a través de las imágenes 2D obtenidas de la sección tisular y el microscopio óptico representa una buena alternativa para recopilar información 3D. Dicha tarea se basa en el emparejamiento (matching) de imágenes 2D de secciones finas empleando la técnica de registro de imágenes [10, 11, 15, 17, 18, 23, 29, 32, 34, 35, 41, 42, 46, 50, 51, 66, 68, 69, 72, 79, 81]. La operación clave para el registro de imágenes en este escenario consiste en compensar la distorsión entre imágenes de secciones consecutivas que se produce tras el corte transversal. Estos cortes son tan delicados como finos (3 a 5 µm). El proceso de preparación (por ejemplo, seccionado, tintado y sellado) puede añadir una variedad de deformaciones no rígidas entre las que se encuentran la flexión, la erosión, la elongación o el desgarro. A resoluciones microscópicas, una mínima deformación es bastante apreciable, y en no pocas ocasiones la precisión es esencial para el éxito del proceso. Con el objetivo de compensar estas deformaciones, el uso de un registro no rígido es esencial, y su éxito depende de la exactitud de las características seleccionadas para emparejar todo el conjunto de imágenes. En aras a fomentar la precisión, la 89 90 Capítulo 5. Registro no rígido de imágenes microscópicas intensidad de los píxeles se encuentra en permanente evaluación, lo que desemboca en un proceso muy costoso computacionalmente, sobre todo si se emplean métricas de comparación ya consolidadas como la Información Mutua (IM). En este capítulo se desarrolla un método de computación de alto rendimiento para generar un emparejamiento preciso de los rasgos o características (features) de las imágenes entre todas las muestras obtenidas a partir del microscopio. Los retos que deben abordarse para seleccionar y correlar el conjunto de imágenes son muy diversos. Destacan entre ellos los siguientes: 1. Entorno rico en características. La calidad de la textura obtenida a través de las imágenes microscópicas representa una oportunidad única para seleccionar una amplia gama de características sobre las que luego cimentar un robusto emparejamiento. La detección tradicional de estas características, por ejemplo, para la localización de las esquinas, genera un abundante número de parámetros que normalmente perjudican el realismo del emparejamiento final. 2. Imágenes a gran escala. Los escáneres actuales de alta resolución pueden generar imágenes con resoluciones de 0.46 µm/píxel (con objetivos de 40x) e incluso superiores. Así, un desplazamiento de 1 mm. supone más de 2,000 píxeles. Por tanto, sin una buena inicialización, el emparejamiento de características se convierte en una tarea inviable. 3. Elevada carga computacional. Un método habitual para el registro de imágenes de alta resolución consiste en aplicar un registro multiescala donde se reduce la resolución para poder procesar la transformación no rígida, y una vez obtenida esta transformación, se aplica las imágenes de mayor resolución. Este proceso no es trivial porque hay que realizar la transformación de las imágenes flotantes en cada escala. Por ejemplo, en este capítulo las imágenes tratadas tienen una resolución de hasta 23 K x 62 K píxeles y, con un factor de escala de dos, la transformación de las imágenes involucra en torno a 12 K x 31 K píxeles como antesala de la etapa final. Para abordar estos retos, se ha desarrollado una alternativa de alta computación para el registro de imágenes, compuesta de un registro rígido rápido para la inicialización, y un algoritmo eficiente y paralelizable para el registro no rígido implementado en procesadores gráficos (GPUs). Dicho enfoque tiene las siguientes ventajas: 1. Implementación en GPU. La idoneidad de la GPU para nuestro propósito se pone de manifiesto en diversos estudios [9, 58] que subrayan su paralelismo masivo, escalabilidad y bajo coste, rasgos todos ellos muy deseables para el 91 registro de imágenes. Adicionalmente, CUDA mejora esta percepción al ofrecer un modelo de programación alternativo que no requiere conocimientos sobre renderizado o gráficos y permite transformar la tecnología de la GPU en un procesador paralelo en cualquier PC. 2. Registro rígido rápido para la inicialización. El algoritmo de registro no rígido se inicializa con un algoritmo de registro rápido. Este algoritmo usa regiones anatómicas destacadas (por ejemplo, vasos sanguíneos) como características de alto nivel, y la transformación rígida se lleva a cabo con un esquema de votado sobre el espacio de transformación euclídea. Este procedimiento es altamente eficiente y preciso para imágenes histológicas, con la ventaja de poderse ajustar tanto las rotaciones como las traslaciones de forma arbitraria. Esto nos permite comenzar desde un gran punto de partida, reduciendo significativamente el espacio de búsqueda del emparejamiento objetivo. 3. Selección de características. Las características para conseguir un emparejamiento preciso en el método de registro no rígido se seleccionan en base a la complejidad de las proximidades, más que en base a su geometría. De esta forma, a la vez que se reduce la carga computacional, se ofrece al usuario una distribución más uniforme de las características. 4. Correlación cruzada normalizada (NCC) 1rápida para un emparejamiento preciso. La búsqueda precisa de características se basa en el método NCC entre los mosaicos vecinos en cada imagen. El cálculo de NCC puede implementarse eficientemente con la transformada rápida de Fourier (FTT), obteniendo ejecuciones más rápidas que métodos como MI. Por otro lado, NCC tiene una interpretación intuitiva que simplifica la selección de los parámetros empleados como umbrales para discriminar entre buenos y malos emparejamientos. 5. Transformación simple de la salida. Para un emparejamiento preciso, los parámetros de la transformación euclídea procedentes de la inicialización rígida se usan para localizar y transformar los vecinos correspondientes, así se evita aplicar la costosa transformación rígida sobre la imagen completa. 6. Paralelización para el emparejamiento. El proceso de emparejamiento preciso es altamente paralelo, prestándose a ser ejecutado en múltiples núcleos, procesadores o clústeres de computación. La implementación en GPU lleva a cabo la parte del algoritmo más exigente computacionalmente: el cálculo de la correlación cruzada para un emparejamiento 1Normalized Cross Correlation. 92 Capítulo 5. Registro no rígido de imágenes microscópicas de imágenes preciso, lo que supone hasta un 60% del tiempo de ejecución total. Los resultados de la ejecución del algoritmo se exponen en sus versiones serie y paralela, tanto en CPU como en GPU, y haciendo un recorrido por diversos parámetros para estudiar la escalabilidad y eficiencia (ver Tabla 5.2). El conjunto de imágenes para las pruebas experimentales (ver Tabla 5.3) procede de dos proyectos de fenotipo cuantitativo: El primer proyecto es un estudio morfométrico basado en el gen retinoblastoma (un supresor de tumores bien conocido) del desarrollo de la placenta de los ratones. En este estudio, se obtienen tres placentas de control y mutantes con el gen Rb eliminado. Cada muestra se ha seccionado a 3 µm y cada sección se ha tintado usando hematoxilina estándar y eosina. Las secciones tintadas se digitalizan usando un escáner de alta resolución Aperio ScanScope con un objetivo de 20x que produce una resolución de 0.46 µm/píxeles. Las seis muestras constituyen más de 3,000 imágenes con las dimensiones típicas de 16 K x 16 K píxeles, ocupando más de tres terabytes de datos sin comprimir. El segundo proyecto es parte de un estudio reciente sobre el microentorno en los tumores del cáncer de mama en ratones. Las imágenes de este estudio son normalmente de 23 K x 62 K píxeles, y alcanzan un tamaño en torno a cuatro gigabytes sin comprimir. El contenido de este capítulo se estructura de la siguiente manera. Desde la Sección 5.1 hasta la Sección 5.3 se resume el método empleado para el registro de imágenes, incluyendo la revisión del algoritmo rápido de registro rígido en la Sección 5.1, una descripción de la selección de características, el emparejamiento preciso de imágenes en la Sección 5.2, y finalmente un repaso de la transformación de la imagen en la Sección 5.3. Las Secciones 5.4 y 5.5 presentan un resumen de la arquitectura GPU y una descripción de nuestra implementación. Los resultados experimentales se muestran en la Sección 5.6, que posteriormente son analizados en la Sección 5.7. 5.1. Registro de imágenes para la inicialización Aunque el desarrollo del algoritmo del registro rígido queda fuera de nuestra jurisdicción, resulta interesante conocer algunos pormenores del mismo para luego asimilar las virtudes de la versión de alto rendimiento que hemos llevado a cabo. La inicialización proporcionada por el registro rígido reduce el área de búsqueda empleada para el emparejamiento preciso en el registro rígido, y este hecho provoca una reducción importante de la carga computacional. 5.1. Registro de imágenes para la inicialización 93 La base del algoritmo de registro rígido reside en el emparejamiento de características de alto nivel, como regiones que corresponden a estructuras anatómicas específicas como los vasos sanguíneos o los conductos mamarios. El empleo de características de alto nivel presenta una multitud de ventajas respecto a los métodos que usan características más primitivas como puede ser la detección de esquinas. En primer lugar, las características de alto nivel son fácilmente segmentables porque corresponden a grandes regiones continuas de píxeles con similitudes cromáticas. Todo lo que se necesita para la extracción es una simple clasificación por color de los píxeles y una secuencia de operaciones morfológicas [71]. Ambas operaciones disponen de implementaciones optimizadas en las librerías comunes de procesamiento de imágenes. Además, a menudo esta extracción puede realizarse con menor muestreo sin perder fidelidad. Por ejemplo, para el cálculo del registro rígido hemos empleado factores de magnificación 5x sobre imágenes 20x. Puesto que el registro rígido únicamente contribuye en la fase de inicialización, donde aún se trabaja con las imágenes magnificadas originales, la pérdida de información por relajar el muestreo no afecta a los resultados finales. En segundo lugar, el número de características de alto nivel se limita habitualmente incluso en imágenes microscópicas para mantener el número de posibles emparejamientos. Finalmente, estas características puede emparejarse a través de parámetros como la forma y el tamaño, que son descriptores globales. Por lo tanto, el proceso no está limitado a búsquedas locales, pudiendo acomodar imágenes con un rango completo de desalineamiento. El uso de características como la forma y el tamaño como criterio para el emparejamiento reduce además la ambigüedad. La cuestión fundamental para cualquier plan de búsqueda de características está en la detección de las discordancias. Con esta finalidad se han desarrollado dos planteamientos. Cualquier par de características de emparejamiento pueden generar una transformación rígida especificada por la rotación del ángulo θy la transformación T si las distancias dentro de la imagen son consistentes. Esto puede ser concebido como gran parte de las transformaciones basadas en los emparejamientos aproximados que deben concentrase alrededor de los parámetros verdaderos en el espacio euclídeo de transformación. Por tanto, la selección de la transformación rígida óptima se puede obtener como el punto de mayor coincidencia dentro de un proceso de votado. La Figura 5.1 muestra un ejemplo de este proceso. En el caso donde pueda reducirse el número total de características presentes en la imagen (por ejemplo, menos de 10), los resultados del proceso de votado pueden ser menos fiable. Por esta razón, se ha desarrollado un enfoque teórico-gráfico basado en el mismo principio para la búsqueda de características [66]. Este algoritmo tiene tiempos de ejecución razonables incluso en implementacio- 94 Capítulo 5. Registro no rígido de imágenes microscópicas (a) (b) (c) (d) (e) (f) Figura 5.1: Registro rígido rápido usando características de alto nivel. nes modestas. Usando Matlab, el parámetro rígido estimado para el ejemplo de la Figura 5.1 se calculó en menos de cuatro segundos en un sistema de características similares a los descritos en la Sección 5.6. 5.2. Extracción y búsqueda de características Una vez descrita la inicialización rígida, toma el relevo el procesamiento del registro no rígido. La corrección de la distorsión no rígida para obtener la precisión necesaria en el fenotipo cuantitativo exige el establecimiento de un gran número de correspondencias precisas. Estas correspondencias deben ser además distribuidas equitativamente a través de las áreas de la imagen que son de interés para obtener una calidad uniforme. Estas consideraciones se remiten a la extracción de características y al emparejamiento, donde los parámetros rígidos de inicialización y el muestreo identifican las regiones correspondientes y son comparadas usando NCC. 5.2.1. Extracción de características La primera cuestión en la extracción de características es la selección de características no ambiguas que resultan candidatas en emparejamientos más precisos y 5.2. Extracción y búsqueda de características 95 específicos. Esto es especialmente importante en el registro no rígido, ya que usar la colección de emparejamientos para interferir en algo sobre la calidad de un emparejamiento individual es difícil debido a la libertad y menor escala de la distorsión no rígida. En este sentido, los emparejamientos en el registro no rígido son de naturaleza local: la única información disponible para juzgar la calidad procede de la vecindad de la característica. En nuestro caso, los mosaicos seleccionados son de un contenido relevante, una mezcla de diferentes tejidos o tejido y fondo que elaboran una apariencia característica y son candidatos a la generación de emparejamientos muy específicos. Y para asegurar dicho procedimiento, los mosaicos seleccionados tienen una varianza que supera un umbral mínimo establecido. De esta forma, la característica de cualquier punto pcon coordenadas [x y]y centrado en la ventana W x W píxeles, debe cumplir 1 W2−1∑ i,j (t(i, j)−t)2≥σ2(5.1) donde tes el patrón, una representación en escala de grises de la ventana de píxeles centrada en pcon valor medio t, y σ2es el umbral de la varianza. Existe la posibilidad de que este umbral pueda ser superado, dando lugar a un resultado ambiguo (considera una plantilla pequeña con la mitad superior blanca y la mitad inferior negra dentro de una plantilla similar pero de mayor tamaño), aunque este tipo de casos apenas se presentan en la práctica. Para mantener una cantidad razonable de características y poder tratar la distribución de características, las características se han muestreado sobre el espacio de las imágenes con un mosaico de W x W. Por ejemplo, en las imágenes de la placenta de 16 K x 16 K, normalmente se crearían mosaicos en el rango de 150-350 píxeles para generar un total de 2025-11236 características posibles, aunque la gran mayoría de las características se descartan por no alcanzar el umbral de la varianza. Con un conjunto de características identificadas en una imagen, sus correspondencias únicas en la imagen siguiente son fáciles de determinar. 5.2.2. Búsqueda de características Para las características de cualquier punto seleccionado p1con coordenadas [x1 y1] en la primera imagen, en primer lugar se selecciona una ventana de píxeles BxB centrada en p1. Esta ventana se transforma a escala de grises y se rota por un ángulo θobtenido de la inicialización del registro rígido. Posteriormente, el área central 96 Capítulo 5. Registro no rígido de imágenes microscópicas de W1xW1píxeles se usa como plantilla p1para identificar p2, lo que sería el punto correspondiente p1en la segunda imagen. El cálculo de Bdepende de θyW1, abarcando suficiente espacio para acomodar el patrón. La coordenada p2se puede estimar usando la siguiente expresión p2′=[x2 y2]=[cos(θ)−sin(θ) sin(θ) cos(θ)][x1 y1]+T(5.2) donde Tcorresponde con el vector de traslación obtenido de la inicialización del registro rígido. Un patrón de W2xW2píxeles centrado en p2′y designado como ventana de búsqueda se toma de la segunda imagen. El NCC se procesa entre la plantilla y la ventana de búsqueda, y el centro del área que se corresponde con el mayor valor NCC se pone en p2. Si este valor máximo sobrepasa un umbral (normalmente 0.8 o superior), el emparejamiento se considera correcto y p1yp2se guardan como correspondencias. La Figura 5.2 refleja este proceso. La selección de W1yW2se basan en la severidad de la deformación y en la capacidad computacional. Empíricamente se ha fijado W2= 2W1. Una deformación superior necesita una tamaño de ventana de búsqueda W2superior. Los resultados experimentales se han evaluado con diferentes tamaños de W1yW2(ver Sección 5.6), con algunos de ellos favorables en CPU y otros en GPU (ver Tabla 5.2). Figura 5.2: Proceso de emparejamiento de características. (a) Imagen inicial. (b) Imagen inicial rotada por el ángulo obtenido en el registro rígido. (c) Región seleccionada dentro de la primera imagen (parche patrón o ventana de características de W1xW1 píxeles). (d) Región de búsqueda dentro de la segunda imagen (ventana de búsqueda de W2xW2píxeles). Entre las medidas de similitud comúnmente utilizadas para la información de la intensidad y que son diferentes a NCC están la raíz cuadrada de la suma de las diferencias (SSD) y MI. SSD no es una buena elección para las imágenes microscópicas porque el contenido tiende a ser discreto (por ejemplo, los límites entre el núcleo celular, el citoplasma y la membrana celular). MI se usa habitualmente como una métrica 5.2. Extracción y búsqueda de características 97 en estrategias de búsqueda por gradiente, pero el histograma conjunto provoca que la búsqueda exhaustiva sea muy costosa computacionalmente. Por tanto, se ha seleccionado NCC porque, además de su robustez en la identificación de similitudes, es altamente eficiente cuando se implementa con la FFT. La robustez es un parámetro clave en la aplicación, ya que resulta muy frecuente encontrar variaciones de intensidad debido a que las secciones no son uniformes en grosor y tintado. Por otro lado, los valores de NCC tienen una interpretación tan intuitiva que la selección de los parámetros umbrales es relativamente fácil. 5.2.3. Procesamiento de NCC La gran cantidad de características que existen dentro de un conjunto de datos convencional hace que el procesamiento eficiente de NCC sea una parte fundamental. Además, más que usar una estrategia de búsqueda, NCC se procesa entre el patrón y las ventanas de búsqueda dentro de todos los desplazamientos posibles con el fin de evitar problemas de mínimos locales. Dado un patrón tde tamaño W1xW1con media t, y ventana de búsqueda scon tamaño W2xW2,W2> W1, el NCC entre tyses el cociente de la covarianza y las varianzas individuales: ρ(u, v) = ∑ x,y {t(x−u, y −v)−t}{s(x, y)−su,v} ({t(x−u, y −u)−t)}2{s(x, y)−su, v}2)1 2 (5.3) donde su,v es la media de la parte de la ventana de búsqueda que se solapa con el patrón en un desplazamiento (u, v). Para calcular los factores de normalización en el denominador se usa el método de suma acumulada [47]. Este método evita los costosos cálculos locales de la media y la varianza de la ventana de búsqueda para la región solapada cuando el patrón se desplaza a través de las posiciones (W1+W2−1)2, reduciendo el número de operaciones de 3W2 2(W1−W2−1)2a aproximadamente 3W2 1. La correlación cruzada no normalizada del numerador se calcula a través del teorema de convolución de la Transforma de Fourier Discreta (DFT,) que relaciona el producto del espectro de la DFT con la convolución circular en el dominio espacial. Respecto a las correlaciones cruzadas, en este capítulo es de especial interés la correlación convencional, por lo que sytse rellenan con ceros hasta el tamaño W1+W2−1 antes de aplicar la transformada para asegurar que en el resultado no existen porciones de solapamiento circular. La Figura 5.3 muestra algunos ejemplos seleccionados aleatoriamente de las regiones emparejadas. 104 Capítulo 5. Registro no rígido de imágenes microscópicas 2. Las glándulas mamarias para estudiar el microentorno del tumor del cáncer de mama [79]. La Tabla 5.3 recoge los detalles de los conjuntos de imágenes empleados. El objetivo en ambos casos es la reconstrucción 3D de los tejidos para su estudio microanatómico. Campo de Área de investigación Zona del Carga Tamaño de Número de estudio y objetivos biomédicos ratón computacional imagen (píxs.) secciones Genética Estudio funcional de un gen Placenta Media 16K x 16K 100 Oncología Tumor de cáncer de mama Mamas Alta 23K x 62K 4 Tabla 5.3: Conjunto de imágenes de entrada empleadas como datos de entrada para el algoritmo de registro. 5.5.2. Plataforma hardware La aplicación que procesa el registro automático ha sido implementada en el sistema descrito en la Sección 1.3.2, en un nodo GPU de visualización donde las características de un procesador de doble núcleo AMD Opteron 2218 han sido combinadas con una GPU de doble zócalo Nvidia Quadro FX 5600 (ver Figura 1.1). La CPU va acompañada de 4 GB de memoria DRAM DDR2 a 667 MHz, mientras que la GPU dispone de 1.5 GB de memoria DRAM GDDR3 a 1600 MHz (ver resto de características en la Tabla 1.2). El sistema en conjunto tiene 7 GB de memoria, un disco duro de 750 GB SATA II a 7200 RPM con 16 MB de memoria caché y una tarjeta InfiniBand para la comunicación exterior. En los experimentos no se ha considerado el tiempo de lectura de los archivos de las imágenes. En cualquier caso, este tiempo se puede ocultar parcialmente solapando las comunicaciones de entrada/salida con el procesamiento en GPU gracias a que la comunicación entre ambos procesadores es asíncrona. 5.5.3. Software La implementación en GPU se ha desarrollado con la versión 1.1 de CUDA, y para aplicar el paralelismo hacia las dos GPUs hemos empleado la librería p-threads (POSIX threads). En la CPU hemos utilizado el compilador de Microsoft Visual Studio 2005 8.0. Además, una implementación previa en Matlab 7.1 ha servido como herramienta para 5.6. Resultados experimentales 105 el prototipado y validación de los resultados, a la vez que como tiempo de ejecución de referencia. 5.6. Resultados experimentales Tal y como refleja la Tabla 5.3, una numerosa cantidad de experimentos se han llevado a cabo sobre 100 imágenes de conjunto de la placenta y sobre cuatro imágenes para el conjunto de las mamas. 5.6.1. Resultados del registro de imágenes Una evaluación directa de la calidad de los algoritmos para el registro de imágenes microscópicas es una tarea complicada debido a dos razones principales: (1) la falta de una referencia ideal y (2) la validación del algoritmo no atañe al contenido de este capítulo. La demostración que se ha realizado sobre el registro de algunas muestras ha sido visual tal y como se muestra en la Figura 5.6. En otro estudio que no hemos presentado aquí se aprecia que la discrepancia entre las imágenes queda reducida a diez píxeles, siendo esta cantidad insignificante en comparación con el tamaño de los datos tratados. Dicha discrepancia está dentro de lo esperado por la diferencia morfológica entre las imágenes. Además de comparar los resultados del registro rígido y no rígido de imágenes, se han analizado las posibilidades para que el método de registro sea útil para el objetivo de la reconstrucción 3D. Como se demuestra en la Figura 5.6f, las imágenes registradas se apilan y se muestran usando una presentación volumétrica. Las secciones cruzadas virtuales generadas demuestran la necesidad del registro no rígido para este caso. Es evidente que los resultados del registro no rígido dirigen la reconstrucción suave de las estructuras microscópicas mientras que el registro rígido dirige los límites sobresalientes de la estructura. 5.6.2. Caracterización de la carga de trabajo Hay que destacar que el tiempo de ejecución para cada imagen dentro del mismo conjunto experimenta variaciones debido al contenido y a que, como consecuencia de ello, el número de características procesadas es diferente. Como se describe en la Sección 5.2.1, la varianza se calcula en una ventana de 200 x 200 píxeles para conservar solamente los puntos que son representativos. Este procedimiento puede provocar que 106 Capítulo 5. Registro no rígido de imágenes microscópicas Figura 5.6: (a) Ejemplo de reconstrucción 3D de la placenta de ratón. Dado que las imágenes son de grandes dimensiones, sólo se muestra una fracción de la reconstrucción correspondiente a 30 secciones. (b-e) Registro de las imágenes de glándulas mamarias: (b) área de 1000 x 1000 píxeles de la imagen de referencias; (c) el correspondiente área de 1000 x 1000 píxeles de la imagen flotante; (d) el área de la imagen flotante después de la transformación no rígida; (e) solapamiento de las dos imágenes. (f) imágenes registradas (arriba las imágenes del registro rígido y abajo las del registro no rígido) apiladas y presentadas usando un renderizado en volumen. La vistas frontales son las secciones cruzadas virtuales generadas tras el apilamiento 3D. imágenes de diferente tamaño generen una carga computacional diferente basada en el contenido (cuanto más homogénea es la imagen, menor es la carga computacional). La Tabla 5.4 resume el número de características extraídas para cada imagen de entrada perteneciente al conjunto de datos mamario, así como el tiempo de ejecución y total para completar el algoritmo de registro en una CPU Opteron. Número de características extraídas Carga de trabajo (en segundos) Tamaño ventana: Pequeño Medio Grande Tiempo de ejeTiempo de eje- (patrón, búsqueda) (342,683) (500,1000) (683,1366) cución con E/S cución sin E/S Mama 1 1196 655 384 650.86 558.54 (85 %) Mama 2 1048 568 312 497.83 414.17 (83 %) Mama 3 3119 1528 854 1320.01 1192.69 (90 %) Mama 4 690 322 168 463.77 340.62 (73 %) Tabla 5.4: Tamaños de las ventanas de búsqueda y características empleadas en los algoritmos de registro (en píxeles). El porcentaje de características procesadas oscila entre 4 % y 30 % del área total de la imagen, variando ligeramente estos valores según el tamaño de la ventana (ver Tabla 5.2). Sin embargo, se puede considerar que los porcentajes son estables para cada imagen si se selecciona el tamaño más pequeño de imagen como la más representativa (resolución de búsqueda mayor). Bajo esta hipótesis, la Figura 5.7 proporciona 5.6. Resultados experimentales 107 los detalles sobre el porcentaje de las características procesadas para el conjunto de imágenes de la placenta y de las mamas. El mínimo porcentaje para las imágenes de la placenta ocurre en la imagen 5 con un 10.48 %, mientras que el máximo se da en la imagen 99 con un 30.38%, y una medida total de 19.88%. En las imágenes de las glándulas mamarias el mínimo porcentaje es 4.82 % para la imagen 4, con un máximo de 20.71% para la imagen 3, y una media de 10.77 %. De acuerdo a nuestra definición de característica, el conjunto de imágenes de la placenta contiene aproximadamente el doble de información relevante, mientras que el conjunto de las glándulas mamarias representa una matriz con una tasa superior de dispersión incluso con un tamaño superior de imágenes. (a) Conjunto de imágenes de la placenta. (b) Conjunto de imágenes mamarias. Figura 5.7: Porcentaje de características procesadas por imagen en cada conjunto de imágenes de entrada. El tamaño para la ventana de búsqueda y de plantilla es el pequeño por ser el más representativo. 5.6.3. Tiempos de ejecución en CPU La Figura 5.8 muestra los tiempos de ejecución para el algoritmo de registro representado en la Figura 5.4 cuando se procesa completamente en la CPU usando la librería FFTW. Los resultados para el conjunto de imágenes de la placenta y las glándulas mamarias se muestran a la izquierda y derecha, respectivamente. Dentro de cada caso, hemos realizado las pruebas con tres patrones y ventanas de búsqueda diferentes (ver Tabla 5.2): pequeñas (azul, más a la izquierda), mediana (rojo, céntrico) y grande (verde, más a la derecha). Según los consejos proporcionados por la librería FFTW, los tamaños pequeño y grande cumplen las condiciones óptimas, mientras que el tamaño mediano rompe todas las reglas. Este hecho ralentiza la ejecución, con un tiempo medio para el caso de la placenta de 294.57 segundos para el tamaño mediano, 57.97 segundos para el tamaño pequeño y 91.33 segundos para el tamaño grande. El resultado supone un incremento del 57% cuando se duplica el tamaño de la ventana dentro de las condiciones óptimas y un 222% cuando no se cumplen dichas condiciones. 108 Capítulo 5. Registro no rígido de imágenes microscópicas El comportamiento para las imágenes mamarias es similar, aunque los incrementos se reducen al 26% y 147 % respectivamente, con tiempos de ejecución medios de 530.41 segundos (tamaño pequeño), 1660.91 segundos (mediano) y 669.96 segundos (grande). (a) Conjunto de imágenes de la placenta. (b) Conjunto de imágenes mamarias. Figura 5.8: Tiempos de ejecución en la CPU Opteron para el algoritmo de registro sobre un par de imágenes de diferentes conjuntos de muestra y tamaños de ventana. El primer par de números de la leyenda corresponde con el tamaño pequeño de ventana de búsqueda y plantilla, respectivamente. Los dos pares siguientes se corresponden con el tamaño mediano y grande, respectivamente. 5.6.4. Tiempos de ejecución en GPU La Figura 5.9 muestra los tiempos de ejecución para el algoritmo de registro cuando la GPU asiste a la CPU procesando la correlación cruzada basada en FFT con CUDA. La gráfica de la izquierda es el resultado para el conjunto de imágenes de la placenta, y a la derecha están los resultados para las muestras mamarias, diferenciándose con leyendas los distintos tamaños de ventana (ver Tabla 5.2). En esta ocasión, los tamaños pequeño y grande cumplen todas las condiciones impuestas por la librería CUFFT, y además el tamaño mediano de 794 píxeles satisface la condición de ser múltiplo de un número primo pequeño (en este caso, 7). Sin embargo, la sobrecarga es aún apreciable. Los tiempos medios para la placenta son de 19.27 segundos (pequeño), 47.80 segundos (mediano) y 22.22 segundos (grande), y la disminución es del 15% cuando el tamaño de ventana asciende al doble dentro de las condiciones óptimas, y un 115% adicional fuera de dichas condiciones. Para las imágenes mamarias, los tamaños grandes consiguen mejorar ligeramente el rendimiento de los pequeños, y la sobrecarga fuera de las condiciones óptimas (tamaño mediano) alcanza el mayor valor: 531 %. 5.6. Resultados experimentales 109 (a) Conjunto de imágenes de la placenta. (b) Conjunto de imágenes mamarias. Figura 5.9: Tiempos de ejecución en la GPU Quadro para el algoritmo de registro sobre un par de imágenes de diferentes conjuntos de muestra y tamaños de ventana. El primer par de números de la leyenda corresponde con el tamaño pequeño de ventana (de búsqueda y plantilla), y los siguientes, con el tamaño mediano y grande. 5.6.5. Comparativa CPU-GPU La fila central de la Tabla 5.5 muestra los factores de aceleración media cuando la GPU procesa la correlación cruzada basada en FFT utilizando CUDA. Las ganancias no son estables para los casos que están fuera de las condiciones óptimas, y los resultados más reales se dan en los tamaños pequeño y grande, donde los tamaños de la ventana siguen estrictamente las indicaciones marcadas por las librerías FFTW y CUFFT. En el caso de la placenta, las imágenes pequeñas generan un factor de aceleración de 3.00x mientras que para las grandes se alcanza 4.11x. En el caso de las glándulas mamarias, dichos factores son más modestos: 2.00x y 2.59x, respectivamente. Imagen de entrada: Placenta: 16K x 16K Mamaria: 23K x 62K Tamaño de ventana: Pequeño Mediano Grande Pequeño Mediano Grande (búsqueda, patrón) (171,342) (250,500) (342,683) (342,693) (500,1000) (683,1366) CPU tiempo ejec. 57.97 294.57 91.33 530.41 1660.91 669.96 GPU tiempo ejec. 19.27 47.80 22.22 264.09 1629.72 257.95 GPU factor acel. 3.00x 6.16x 4.11x 2.00x 1.01x 2.59x Tiempo 2 GPUs 13.13 26.05 13.66 225.17 837.51 234.62 2 GPU / 1 GPU 1.46x 1.83x 1.62x 1.17x 1.94x 1.09x 2 GPU / 1 CPU 4.41x 11.30x 6.68x 2.57x 1.98x 2.85x Tabla 5.5: Tiempos de ejecución (en segundos) y factores de aceleración para las diferentes implementaciones desarrolladas para el algoritmo de registro con máximo rendimiento. La Figura 5.10 demuestra que el factor de mejora en la GPU tiene gran dependencia de la imagen de entrada, especialmente para el conjunto de imágenes mamarias 110 Capítulo 5. Registro no rígido de imágenes microscópicas donde los números son menos consistentes. Adicionalmente, estas ganancias son más inestables conforme aumenta el tamaño de las ventanas. Este efecto se debe a que el contenido de las imágenes más grandes empieza a ser más heterogéneo en grandes búsquedas, mostrando por otro lado disparidades entre las imágenes. La Figura 5.8 corrobora dicho efecto. (a) Conjunto de imágenes de la placenta. (b) Conjunto de imágenes mamarias. Figura 5.10: Comparativa entre los tiempos de ejecución de la CPU y la GPU en términos de factores de aceleración. Cuando el tamaño de las ventanas aumenta, los tiempos de ejecución son más irregulares en (b). El primer par de números de la leyenda corresponde con el tamaño pequeño de ventana (de búsqueda y plantilla), y los siguientes, con el tamaño mediano y grande. (a) El factor medio para la placenta es de 3.00x (pequeño), 6.16x (mediano) y 4.11x (grande). (b) Para las mamas el factor medio es de 2.00x (pequeño), 1.01x (mediano) y 2.59x (grande). 5.6.6. Paralelismo y escalabilidad en la GPU La popularidad conseguida por la GPU durante la última década se debe a la espectacular escalabilidad de su arquitectura, manteniéndose en su objetivo de duplicar su rendimiento cada seis meses. Además de la tendencia intra-chip, otras iniciativas como la tecnología SLI de Nvidia y Crossfire de ATI han surgido para explotar el paralelismo inter-chip (SMP - MultiProcesamiento Simétrico). La iniciativa ha alcanzado un notable éxito dentro de la industria de los videojuegos, pero su impacto sobre la comunidad GPGPU ha sido bastante menor. Esta sección evalúa el rendimiento del algoritmo de registro sobre un par de GPUs cuando se aplica paralelismo SMP. Nuestras técnicas de programación son extensibles a cualquier número de tarjetas gráficas, ya que los métodos usados para particionar los datos garantizan perfectamente la escalabilidad del problema. Sin embargo, en este ambicioso algoritmo hay que advertir el papel crítico que asume el sistema de entrada/salida: docenas o incluso cientos de GPUs que trabajan en paralelo pueden encontrar un método fácil para distribuir las diferentes ventanas de búsqueda de una 5.6. Resultados experimentales 111 manera eficiente cuando tratan con grandes imágenes, pero el sistema de archivos tendría que ser de tal rendimiento que permitiese leer cada imagen en paralelo con un ancho de banda sostenible y suficiente para conseguir una tasa de procesamiento cercana al Teraflop. Durante los experimentos, este cuello de botella no ha sido analizado en detalle, pero la Tabla 5.4 cuantifica en las dos últimas columnas el tiempo de ejecución total (incluyendo entrada/salida) y el tiempo de cómputo (excluyendo entrada/salida), revelando así que la interacción con disco es responsable del 10-20% del tiempo total de ejecución. Este tiempo no se incluye en los análisis posteriores, ya que es el mismo tanto para la CPU como para la diferentes versiones de la implementación en GPU y no es el objeto de este trabajo adentrarnos en el sistema de almacenamiento. Implícitamente se asume que los datos de las imágenes se encuentran en la memoria DRAM o que se pueden recuperar de una manera eficiente utilizando un sistema de archivos en paralelo o en RAID. Una vez que los datos se encuentran en la CPU, hay dos alternativas básicas para distribuir la carga computacional del algoritmo de registro entre múltiples GPUs: por bloques o cíclica. Para el caso concreto de un par de GPUs, pero sin perder generalidad, la distribución por bloques asigna la mitad superior de la imagen a una GPU y la mitad inferior a la otra. La distribución cíclica es justo al contrario, numera las diferentes imágenes para enviar las impares a una CPU y las pares a la otra. Dado que las características más interesantes de las imágenes tienden a estar concentradas, la distribución por bloques es más congruente y ha sido la seleccionada para nuestros experimentos. Nuestro método de paralelización funciona de la siguiente forma: se crea un hilo para cada región de la imagen que calcula la varianza en la CPU en cuestión para evaluar si merece la pena procesarla. Si se supera esta prueba, la imagen se envía a una GPU predeterminada para procesar la correlación cruzada normalizada y la búsqueda de características. La Tabla 5.6 resume el número de imágenes procesadas y descartadas en cada GPU dependiendo de la imagen de entrada procedente del conjunto de imágenes mamarias. El desequilibrio de la carga computacional está entre un 2.76% de la imagen 2 y un 13.33% de la imagen 4, en sentido creciente conforme es menor el número de imágenes a procesar (tasa de dispersión de la imagen de entrada). Finalmente, la Figura 5.11 desvela que la ganancia obtenida cuando una segunda GPU entra en juego es muy diversa, empezando con 30-50% para tamaños pequeños de ventana, continuando con 60% para tamaños grandes de ventana, y acabando con una escalabilidad óptima (100 %) en tamaños medianos. Estas ganancias son proporcionales a la carga computacional, demostrando que la GPU es más escalable cuanto más puede explotar su intensidad aritmética. O lo que es lo mismo, los GFLOPS no 112 Capítulo 5. Registro no rígido de imágenes microscópicas Imagen de ProceNúmero de mosaicos Desequilibrio de la Tiempo de entrada sador Testeados Procesados/descartados carga de trabajo ( %) ejec. (s.) Mama 1 GPU 1 1672 196/1476 4.08 260.41 GPU 2 1672 188/1484 Mama 2 GPU 1 1496 158/1338 2.53 101.32 GPU 2 1496 154/1342 Mama 3 GPU 1 1872 428/1444 2.76 522.43 GPU 2 1911 426/1485 Mama 4 GPU 1 1786 78/1708 13.33 225.37 GPU 2 1786 90/1696 Tabla 5.6: Número de mosaicos procesados y descartados para cada imagen en el conjunto de imágenes mamarias en cada GPU bajo la ejecución paralela en dos GPUs. están limitados por la escasez de datos procedente de un ancho de banda insuficiente entre la memoria de vídeo y la GPU. (a) Conjunto de imágenes de la placenta. (b) Conjunto de imágenes mamarias. Figura 5.11: Escalabilidad de la GPU: factores de mejora cuando se habilita una segunda GPU. El primer par de números de la leyenda corresponde con el tamaño pequeño de ventana (de búsqueda y plantilla), y los siguientes, con el tamaño mediano y grande. (a) El factor medio para la placenta es de 1.46x (pequeño), 1.83x (mediano) y 1.62x (grande). (b) Para las mamas el factor medio es de 1.17x (pequeño), 1.94x (mediano) y 1.09x (grande). 5.7. Conclusiones Tras un análisis de los experimentos se pueden extraer la siguientes conclusiones: 1. El conjunto de imágenes de la placenta muestra factores de aceleración superiores en la plataforma gráfica gracias a que las imágenes tienen un contenido más relevante, conducente a una mayor carga computacional que explota mejor 5.7. Conclusiones 113 la intensidad aritmética y el ancho de banda. Además, un numero pequeño de características procesadas implica una alta presencia de sentencias condicionales en el código, uno de los rasgos más perjudiciales para el rendimiento de la GPU. 2. La mayor escalabilidad se consigue con el conjunto de imágenes de la placenta, y sus ganancias son más estables entre los diferentes tamaños de ventana. La gran disparidad de las imágenes mamarias juegan un papel negativo en la distribución de la carga, introduciendo desbalances e impidiendo que el paralelismo sea completamente explotado. En conjunto, la GPU consigue un factor de aceleración de 3-4x en la mayoría de los escenarios típicos (texto enmarcado en Tabla 5.5) comparado con la CPU, y un par de GPUs muestra una escalabilidad satisfactoria, aunque ganancias inestables bajo diferentes imágenes y tamaños de ventana. 120 Capítulo 6. Optimizando los momentos de Zernike sobre Kepler seno) y su distancia, distinguiendo además los píxeles que quedan dentro y fuera del círculo unidad. ésta es la fase que aprovecha la optimización de simetría. 2. Polinomios de Zernike. Con las distancias y senos/cosenos procedentes de la fase anterior, se calculan los polinomios de Zernike para cada píxel. El sumatorio de la ecuación 6.5 se realiza sobre un espacio de memoria compartida para evitar accesos reiterados a memoria global de CUDA. 3. Aplicación a la imagen de entrada. El resultado del espacio conseguido en la fase anterior se multiplica con la imagen de entrada. 4. Sumatorio de píxeles. Se contabiliza la suma de los píxeles comprendidos dentro de círculo unitario a través de un algoritmo de reducción. 5. Sumatorio de las componentes de cada píxel. Al igual que en la fase previa, sumamos la componente que cada píxel aporta al momento de Zernike, fusionándola en un único valor mediante un algoritmo de reducción. Con esta implementación de partida en CUDA, los experimentos son realizados el servidor Yuca descrito en la Sección 1.3.3, la Tabla 1.3 presenta el hardware sobre el que realizaremos todos los experimentos. Contaremos con dos GPUs de generaciones diferentes: Una Fermi GF100 que sitúa los tiempos de referencia para la arquitectura antecesora, y otra Kepler GK110 que permite cuantificar las mejoras logradas en los nuevos procesadores SMX con paralelismo dinámico y Hyper-Q. 6.4. Optimizando Zernike sobre Kepler En esta sección analizaremos las distintas partes del código en las que se pueden aplicar las características de la arquitectura Kepler, aunque como veremos no todas ellas revertirán en mejoras productivas. 6.4.1. Recursividad Comenzamos describiendo la implementación del método recursivo en GPU, a pesar de que pueda parecer más desafiante que el directo para lograr una ejecución eficiente [48]. Dentro de los métodos recursivos, q-recursive es el más actual y eficiente [16] para computar todas las repeticiones de los momentos de Zernike que corresponden a un orden concreto. Los dos primeros momentos calculados corresponden a las dos repeticiones más altas, y a partir de ahí, se obtienen todas las repeticiones progresivamente 6.4. Optimizando Zernike sobre Kepler 121 más bajas a partir de unas expresiones estáticas que involucran a los dos momentos de repetición inmediatamente superiores. Esta metodología resulta acertada en CPU cuando el objetivo es computar varios momentos pertenecientes a un mismo orden, aunque para ello es necesario cambiar algunos aspectos de la implementación base que comentamos en la Sección 6.3. La codificación de este método se va a llevar a cabo siguiendo un proceso iterativo que calcula en cada paso los polinomios de Zernike para una repetición dada a partir de los momentos previamente almacenados. El espacio de memoria aumentará en el mismo factor que el número de repeticiones a calcular, pero la complejidad del algoritmo disminuye y los kernels que no dependen de la repetición se pueden amortizar para todas las iteraciones. Los kernels numerados como 1 y 4 en la Sección 6.3 permanecen intactos, mientras los demás sufren los cambios que se detallan a continuación: El kernel 2 que aplica los polinomios de Zernike a cada píxel debe escindirse en dos, ya que la naturaleza del algoritmo iterativo impide enlazar el procesamiento de los polinomios de Zernike con su aplicación al punto del espacio de coordenadas cartesianas. Ahora tenemos: 2.1 Polinomios de Zernike en forma recursiva. Procesa los polinomios de Zernike recursivamente. Este kernel se ejecuta tantas veces como repeticiones haya para los momentos de Zernike de un orden específico. Cada valor procesado usa su propio espacio de memoria. 2.2 Aplicación al espacio cartesiano. Los polinomios de Zernike se aplican sobre todo el espacio, junto con las funciones trigonométricas que le anteceden. Los kernels 2.1, 3 y 5 aumentan su carga de trabajo en el mismo factor que el número de repeticiones, recayendo este trabajo sobre cada hilo, que distingue de forma unívoca la partición del espacio de memoria al que necesita acceder gracias a su identificador de bloque e hilo dentro de éste. 6.4.2. Paralelismo dinámico El paralelismo dinámico puede aplicarse de varias formas a los momentos de Zernike según la carga de trabajo que se traslade desde la CPU a la GPU. A continuación se exponen diferentes estrategias de paralelismo dinámico que se podrían combinar, y lo que cada una de ellas puede aportar: 122 Capítulo 6. Optimizando los momentos de Zernike sobre Kepler 1. Lanzar kernels desde la GPU. Las llamadas de los cinco kernels del método directo se trasladan al ámbito de la GPU, de forma que inicialmente se lanza un kernel de un único hilo y bloque. Este kernel inicial o raíz realiza las llamadas correspondientes al código de los momentos de Zernike desde la GPU al igual que lo hacía la versión convencional desde la CPU. 2. Calcular una repetición desde cada hilo. En este caso se aprovecha el paralelismo dinámico para calcular todos los momentos de un orden dado. En la versión convencional y siguiendo el método directo, sería necesario aplicar un bucle que itere tantas veces como repeticiones existan. Aplicando paralelismo dinámico, la implementación sería equivalente a la del punto anterior, con la salvedad de que el kernel no tendría un solo hilo, sino uno por cada repetición. Para evitar redundancia en los cálculos, los kernels en común para todas las repeticiones serían procesados desde un mismo hilo. 3. Paralelizar el bucle de los polinomios de Zernike. El bucle for que se requiere en el cálculo de los polinomios de Zernike por el método directo (ver Figura 6.1), es buen candidato para aprovechar el paralelismo dinámico. Cada píxel precisa de dicho cálculo, así que cada uno de ellos lanzará un nuevo kernel que procese de forma concurrente el trabajo de ese bucle. 6.4.3. Hyper-Q El aprovechamiento de Hyper-Q requiere establecer un proceso con varios flujos de ejecución concurrente que no presenten dependencias. En los momentos de Zernike se consigue, al igual que el segundo punto del anterior apartado, cuando necesitamos calcular todas las repeticiones de los momentos de Zernike para un orden específico. El aprovechamiento de Hyper-Q es transparente al programador. La implementación se lleva a cabo usando streams que puedan beneficiarse de las múltiples colas de ejecución concurrente. El aumento del número de colas de 16 a 32 resulta también clave en la nueva arquitectura para aprovechar al máximo el creciente número de procesadores CUDA, ya que ahora, con un número cercano a los tres mil, resulta más probable que un kernel sólo pueda ocupar una fracción de éstos. Distribuyendo los kernels en streams siempre que sea posible, conseguiremos que el remanente de procesadores que haya dejado libre la malla de hilos definida para un primer kernel en ejecución, pueda ser aprovechada desde otros procedentes de streams adicionales. 6.5. Resultados experimentales 123 6.5. Resultados experimentales La Tabla 1.3 resume las prestaciones de las GPUs utilizadas durante la evaluación experimental de nuestras optimizaciones. Hemos empleado un tipo de datos de simple precisión y tomado tamaños de imagen progresivos desde 64x64 hasta 2Kx2K píxeles. El orden máximo de los momentos que se calcula es de 34 debido al límite impuesto por el hardware para el cálculo de los factoriales que aparecen en los polinomios de Zernike. Momentos GPU Tamaño de imagen Mejora de Zernike 64 128 256 512 1024 2048 Mín. Máx. A4,∗ Fermi 0,12 0,17 0,37 1,11 4,08 15,75 0,67x 2,20x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A8,∗ Fermi 0,20 0,32 0,74 2,36 8,83 34,71 0,71x 2,35x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A12,∗ Fermi 0,29 0,50 1,21 4,08 15,32 60,41 0,71x 2,45x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A16,∗ Fermi 0,38 0,72 1,82 6,14 23,49 92,05 0,72x 2,49x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A20,∗ Fermi 0,51 0,97 2,50 8,68 33,37 130,93 0,76x 2,54x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A24,∗ Fermi 0,61 1,27 3,31 11,76 45,39 176,51 0,76x 2,57x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A28,∗ Fermi 0,74 1,59 4,23 15,20 58,20 229,20 0,76x 2,61x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A32,∗ Fermi 0,87 2,00 5,27 18,89 73,30 288,76 0,76x 2,64x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 A34,∗ Fermi 0,95 2,16 5,89 20,96 81,29 319,76 0,78x 2,64x Kepler 0,12 0,17 0,37 1,11 4,08 15,75 Tabla 6.1: Tiempos de ejecución (en milisegundos) para procesar todas las repeticiones de un orden a través del método directo de los momentos de Zernike. Las imágenes de entrada son cuadradas y de dimensiones potencia de dos en un rango entre 64 y 2048 píxeles. 6.5.1. Cambio de arquitectura: SMX En primer lugar, cuantificamos las mejoras producidas por el cambio de arquitectura sin realizar modificaciones en el código. La Tabla 6.1 recoge los tiempos de ejecución sobre Fermi y Kepler para el algoritmo de partida de la Sección 6.3, a la vez que se comparan con el factor de ganancia. El factor de mejora para el cambio de plataforma oscila entre 0.67x y 2.64x. La 124 Capítulo 6. Optimizando los momentos de Zernike sobre Kepler cota mínima de 0.67x, que supone un aumento de tiempo de ejecución de 1.49x, se produce para el tamaño más pequeño de la imagen, y de igual forma, la aceleración máxima de 2.64x llega cuando la carga de trabajo alcanza las dosis más elevadas. La GPU como procesador, y su modelo de paralelismo de datos como forma de programación, se lucen más a medida que crecen las imágenes de entrada, lo que explica que la migración a Kepler acelere más conforme aumenta el número de píxeles a procesar. El tamaño de imagen más pequeño consta de 64 ×64 = 4096 píxeles que se distribuyen en bloques repartidos equitativamente entre todos los multiprocesadores disponibles. El multiprocesador SM de Fermi puede procesar hasta 2 warps de forma concurrente, lo que da un total de 896 píxeles a procesar en nuestra plataforma (2x32x14). El SMX de Kepler llega hasta 4 warps para un total de 1664 píxeles a la vez para aprovechar los recursos. Con estos datos, la imagen pequeña de 64 x 64 píxeles aprovecharía sólo el 32.6% y 18.9% de los recursos para Fermi y Kepler, respectivamente. En este caso, aunque la arquitectura haya sido mejorada, sus recursos son infrautilizados y es la frecuencia, muy superior en Fermi, la que resulta más determinante en el tiempo de ejecución. Con un tamaño de 128 x 128 píxeles, la carga de trabajo ya satura a la máxima que puede procesar concurrentemente Fermi, mientras que la GPU Kepler se queda en un 75.7%. En este tamaño de imagen, el rendimiento tiende a igualarse en ambas arquitecturas, ya que aunque la GPU Kepler no aproveche el 100 % de sus recursos, puede planificar todo el cálculo de una sola vez, mientras que la GPU Fermi necesita una segunda tanda para planificar el remanente de bloques que están a la espera. A partir de aquí, para imágenes de tamaño superior, la mejora de rendimiento en Kepler aumenta hasta alcanzar su cota máxima con el 100 % de los recursos ocupados en esta nueva arquitectura. En un análisis vertical, el aumento del orden en los momentos de Zernike supone una ejecución adicional del algoritmo por cada dos órdenes. Esta carga adicional aumenta el tiempo de cómputo a la vez que beneficia ligeramente a la GPU. 6.5.2. Configuración de la carga de trabajo El incremento del número de núcleos del multiprocesador SMX en un factor de seis respecto a SM predispone a reflexionar sobre si la configuración óptima del tamaño del bloque de hilos puede haber variado respecto a la implementación base. En la GPU Fermi dotada de 32 núcleos, un valor entre 128 y 256 hilos era el óptimo para lograr un factor de ocupación del 100% mientras no hubiera restricciones respecto al uso de los registros y la memoria compartida. La GPU Kepler, cuyo multiprocesador SMX cuenta con 192 núcleos, merece un análisis más pormenorizado a este respecto. 6.5. Resultados experimentales 125 Figura 6.2: Tiempo de ejecución en Kepler para evaluar el rendimiento (escala derecha) en función del tamaño de bloque. Estos valores se contrastan con los teóricos de ocupación para la arquitectura (escala izquierda). La Figura 6.2 desvela los tiempos de ejecución para el kernel que procesa los polinomios de Zernike, siendo éste el más critico en cuanto a recursos. Los tiempos, según la escala de la derecha, se contrastan frente al valor teórico para bloques con diferentes configuraciones de hilos, lo que desemboca en el factor de ocupación según la escala de la izquierda. Aunque teóricamente el mejor rendimiento se consigue para potencias de dos a partir de 128 hilos por bloque, el resultado ha sido levemente mejor para una configuración exacta de 128 hilos, la más baja entre las candidatas. Sin embargo, la mayor divergencia entre la teoría y la práctica se produce cuando el número de hilos por bloque no es suficiente para aprovechar los recursos, o también cuando al añadir un nuevo bloque se sobrepasa el límite de hilos por multiprocesador y, por tanto, hay que sacrificar un bloque en la ejecución concurrente. Ambas zonas se pueden diferenciar: 1. Si el bloque tiene menos de 128 hilos, la limitación para no llegar al número de hilos máximo por multiprocesador la impone el máximo número de bloques permitido. En la GPU Kepler, se pueden ejecutar hasta 16 bloques activos en cada multiprocesador SMX, y la fórmula de la ocupación tiene la siguiente expresión: Ocupacio′ n=Hilos ∗16 2048 (6.8) 2. Cuando el número de hilos alojados en un multiprocesador se acerca a su valor máximo, la inclusión de un nuevo bloque puede sobrepasar el máximo y, por 126 Capítulo 6. Optimizando los momentos de Zernike sobre Kepler tanto, el número permitido de bloques no consigue aprovechar todos los recursos. El peor caso (69 % de ocupación) se da con una configuración de 704 hilos por bloque cuando se aplica la siguiente expresión: Ocupacio′ n=Hilos ∗Bloquesmax 2048 (6.9) 6.5.3. Paralelismo dinámico Las estrategias de paralelismo dinámico descritas en la Sección 6.4.2 van a permitir conocer cuándo es beneficiosa su aplicación. En la Tabla 6.2 se muestra que la aplicación más sencilla, denominada “Lanzar kernels desde la GPU”, y la estrategia “Calcular una repetición desde cada hilo” no resultan positivas. Estos resultados se han tomado para momentos con orden y repetición específicos para, además de comparar los resultados con paralelismo dinámico, dar una idea sobre los tiempos de ejecución para momentos individuales. En el primer caso, los tiempos de ejecución son mayores en un factor entre 1.15x y 2.15x (ver Tabla 6.2a). Esta varianza depende de la carga de trabajo que supone cada llamada a un nuevo kernel desde la GPU. Las imágenes de tamaño superior sufren menos penalización al igual que los momentos computacionalmente más exigentes, que son aquellos en los que la diferencia entre el orden y la repetición es mayor. Para la segunda estrategia, el rendimiento empeora hasta en un factor de 1.4x en el mismo sentido que la estrategia anterior (ver Tabla 6.2b). En este caso, el rendimiento es mejor para las imágenes pequeñas debido a que se procesa un conjunto de momentos que consiguen aprovechar los recursos, mientras que para momentos específicos no sería posible. En resumen, el rendimiento para imágenes pequeñas es sensible a dos aspectos: el propio de lanzar los kernels internamente desde GPU, y el que permite paralelizar los cálculos de los momentos que forman parte de la cadena An,∗. La tercera estrategia, denominada “Paralelizar el bucle de los polinomios de Zernike”, es muy pretenciosa por su naturaleza dinámica. Sin embargo, las pruebas experimentales dejan ver rápidamente que su rendimiento es catastrófico. La implementación requiere que cada hilo lance un nuevo kernel, lo que dispara el tiempo hasta factores 1000x rápidamente. Hemos medido el tiempo que consume una nueva llamada a GPU, resultando entre 5 y 16 µsg., mientras que en la CPU consume alrededor de 3 µsg. No parece lógico asumir que el lanzamiento de un kernel introduzca mayor sobrecarga cuando se lanza desde circuitería mucho más cercana como la propia GPU, y pensamos que esta rémora se solventará en versiones más maduras de los drivers 6.5. Resultados experimentales 127 Momento Sin P. Dinámico Con P. Dinámico Ganancia de Zernike 64 256 1024 64 256 1024 Min. Max. A0,00,08 0,10 0,56 0,16 0,20 0,88 0,48x 0,64x A6,20,08 0,13 0,93 0,16 0,22 1,26 0,53x 0,74x A12,00,09 0,16 1,43 0,16 0,27 1,79 0,58x 0,80x A25,13 0,09 0,17 1,52 0,16 0,27 1,88 0,59x 0,81x A34,00,12 0,28 3,07 0,18 0,38 3,51 0,65x 0,88x A34,18 0,10 0,19 1,82 0,16 0,29 2,20 0,60x 0,83x A34,34 0,08 0,10 0,63 0,16 0,20 0,95 0,49x 0,66x (a) Tiempos de ejecución en milisegundos sin/con paralelismo dinámico. Todos los momentos Tamaño de imagen para un orden dado 64 128 256 512 1024 2048 A4,∗0,73 0,70 0,78 0,70 0,71 0,72 A8,∗0,82 0,76 0,71 0,75 0,76 0,76 A12,∗0,89 0,82 0,76 0,79 0,79 0,79 A16,∗0,91 0,84 0,78 0,80 0,81 0,81 A20,∗0,93 0,86 0,81 0,82 0,82 0,82 A24,∗0,93 0,87 0,82 0,84 0,84 0,84 A28,∗0,95 0,89 0,84 0,85 0,85 0,85 A32,∗0,96 0,91 0,85 0,86 0,86 0,86 A34,∗0,96 0,90 0,85 0,86 0,86 0,86 (b) Factores de aceleración logrados con paralelismo dinámico. Tabla 6.2: Tiempos de ejecución y factores de aceleración logrados para las diferentes estrategias que explotan el paralelismo dinámico con imágenes cuadradas de diferentes tamaños para momentos de Zernike concretos en el primer caso, y todas las repeticiones de cada orden en el segundo. y/o implementaciones más maduras de los multiprocesadores SMX. No obstante, el paralelismo dinámico está orientado a aplicaciones con una filosofía “divide y vencerás” y debe concebirse en la misma línea que la computación con GPUs: Potenciando pocas llamadas a kernels con gran cantidad de datos. Otra característica, que limita el ámbito de aplicación del paralelismo dinámico, es la restricción de que cada hilo no puede acceder a la memoria compartida de los kernels padres. La información a compartir entre kernels padres e hijos queda relegada al uso de la memoria global, con la penalización de rendimiento que esto supone. 6.5.4. Hyper-Q Para aprovechar Hyper-Q en nuestro algoritmo definiremos un stream por cada repetición cuando se persigue calcular todas las repeticiones de un orden. La Figura 128 Capítulo 6. Optimizando los momentos de Zernike sobre Kepler (a) Mejoras logradas con Hyper-Q en Kepler. (b) Mejoras logradas con Kernels concurrentes en Fermi. Figura 6.3: Ganancia obtenida con el uso de streams para distintos tamaños de imagen cuando aumentamos el conjunto de momentos de Zernike a procesar. 6.3a compara el factor de ganancia cuando entra en escena Hyper-Q para este supuesto, variando el orden de los momentos de Zernike y el tamaño de la imagen de entrada sobre la que se aplica. Mientras que la ganancia máxima obtenida es de 2.2x, la ganancia mínima queda a la par con la implementación básica. Esta situación de paridad se produce ya para una carga de trabajo superior a la que abastece a la totalidad de recursos de la arquitectura. Si la imagen es de tal tamaño que la carga de trabajo generada mantiene a la GPU completamente ocupada, no quedan remanentes para ser aprovechados desde streams adicionales mediante Hyper-Q, y cada nuevo stream acabará serializándose en la cola de trabajo. Por otro lado, en los tamaños de imagen más pequeños, la ganancia aumenta proporcionalmente al número de repeticiones a calcular por cada momento. Este escenario es el ideal para el aprovechamiento de Hyper-Q, ya que la poca carga de trabajo suministrada para cada imagen se compensa dando entrada al procesamiento de otras imágenes en paralelo. Por tanto, el rendimiento máximo se consigue procesando la imagen más pequeña y el momento de Zernike para el orden más alto. Como el uso de Hyper-Q no supone ningún cambio para el desarrollador y la gestión de colas para procesar los streams es transparente, la misma aplicación ejecutada en la arquitectura Fermi nos proporciona una visión de la mejora que supuso la inclusión de Kernels Concurrentes. La Figura 6.3b muestra, al igual que se hizo en Kepler, el factor de ganancia obtenido por la incorporación de esta técnica en Fermi. Los valores son similares cualitativamente con la peculiaridad de que, para el tamaño 6.5. Resultados experimentales 129 de imagen inferior donde la ganancia es mayor, la GPU llega a su plena ocupación con órdenes de momento inferiores. Concretamente, el valor máximo 1.86x se consigue para 11 streams, situación que corresponde para el cálculo del momento de orden 20 (A20). Hyper-Q ha supuesto en los momentos de Zernike una aceleración máxima de 1.74x respecto a lo ya conseguido con Kernels Concurrentes en Fermi. El análisis anterior engloba conjuntamente la aplicación de Hyper-Q y el uso pleno de los recursos de los multiprocesadores. Para aislar la aceleración correspondiente a Hyper-Q, hemos realizado un experimento con una imagen de 16 x 16 píxeles (que corresponde con nuestro tamaño del bloque de hilos en CUDA). De esta forma, la ejecución convencional de todas las repeticiones de un mismo orden se realiza secuencialmente y sin aprovechar todos los recursos al enviar un único bloque por iteración. Cuando se usa Hyper-Q o Kernels concurrentes, el uso de los recursos se incrementa por la ejecución de las distintas repeticiones en paralelo. Figura 6.4: Análisis comparativo del beneficio conseguido en la GPU cuando la imagen es de 16 x 16 píxeles, correspondiente a un bloque de hilos CUDA, con objeto de aislar la aceleración atribuida a Hyper-Q y Kernels Concurrentes. La Figura 6.4 muestra los resultados de la comparativa para Fermi y Kepler. Los tiempos cuando no se habilita la paralelización con streams favorecen a la GPU Fermi al tratarse de una imagen de entrada muy pequeña, tal y como ya nos ocurrió y explicamos en la Sección 6.5.1 para imágenes de 64x64 píxeles. Cuando se habilita