Spectral-spatial classification of n-dimensional images in real-time based on segmentation and mathematical morphology on GPUs

UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Departamento de Electrónica y Computación Centro de Investigación en Tecnoloxías da Información (CiTIUS) Tesis Doctoral Spectral-spatial classification of n-dimensional images in real-time based on segmentation and mathematical morphology on GPUs Presentada por: Pablo Quesada Barriuso Dirigida por: Dra. Dora Blanco Heras Dr. Francisco Argüello Pedreira Julio de 2015 Dra. Dora Blanco Heras, Profesora Titular de Universidad del Área de Arquitectura de Computadores de la Universidad de Santiago de Compostela Dr. Francisco Argüello Pedreira, Profesor Titular de Universidad del Área de Arquitectura de Computadores de la Universidad de Santiago de Compostela HACEN CONSTAR: Que la memoria titulada Spectral-spatial classification of n-dimensional images in real-time based on segmentation and mathematical morphology on GPUs ha sido realizada por D. Pablo Quesada Barriuso bajo nuestra dirección en el Departamento de Electrónica y Computación y en el Centro Singular de Investigación en Tecnoloxías da Información de la Universidad de Santiago de Compostela, y constituye la Tesis que presenta para obtar al título de Doctor. Y autorizan la presentación de la tesis indicada, considerando que reune los requisitos exigidos en el artículo 34 de la regulación de Estudios de Doctorado, y que como directores de la misma no incurren en las causas de abstención establecidas en la ley 30/1992. Santiago de Compostela, Julio de 2015 Dora Blanco Heras Directora de la tesis Francisco Argüello Pedreira Director de la tesis To I. L. V. My star, my perfect silence. “If you spend too much time thinking about a thing, you’ll never get it done. Make at least one definitive move daily toward your goal.“ Bruce Lee. “With great power there must also come – great responsibility!.“ Amazing Fantasy #15 – The first Spider-Man story. Acknowledgments I would like to express my gratitude to all the people who have supported me during my thesis. Special thanks goes to my thesis advisors, Dora Blanco Heras and Francisco Argüello Pedreira, for their guidance, support and patience. Their continuous confidence in my work and their motivation have been decisive in the completion of this dissertation. I would also like to acknowledge the Department of Electronics and Computer Science, specially to the Computer Architecture Group, and to the Centro Singular de Investigación en Tecnoloxías da Información (CiTIUS) at the University of Santiago de Compostela, for providing the necessary resources and technical support required by this project. My gratitude also goes to Prof. Jon Atli Benediktsson, University of Iceland, and Prof. Lorenzo Bruzzone, University of Trento, for their advise and kind support during my stage in their research groups. My deepest thanks to my family and close friends for their never ending patience and encouragement during hard moments through all these years. Finally, I am also thankful to the following institutions for the funding provided for this work: Ministry of Science and Innovation, Government of Spain, cofounded by the FEDER funds of European Union, under contract TIN 2010-1754, and by Xunta de Galicia under contracts 08TIC001206PR and 2010/28. Julio de 2015 xvi List of Acronyms DWT discrete wavelet transform EAP Extended Attribute Profile EM Expectation Maximization EMAP Extended Multi-Attribute Profile EMP Extended Morphological Profile FE Feature Extraction FPGA Field Programmable Gate Array FS Feature Selection GPU Graphics Processing Unit HPC High Performance Computing HSEG Hierarchical Image Segmentation ICA Independent Component Analysis MM Mathematical Morphology MNF Minimum Noise Fraction MP Morphological Profile MV Majority Vote nDn-dimensional NPP NVIDIA Performance Primitives NWFE Non-parametric Weighted Feature Extraction OA Overall Accuracy OAO One-Against-One PC Principal Component xvii PCA Principal Component Analysis PR post-regularization RBF Radial Basis Function RCMG Robust Color Morphological Gradient ROSIS-03 Reflective Optics System Imaging Spectrometer SE structuring element SM Streaming Multiprocessor SP Scalar Processor SVM Support Vector Machine WT Wavelet Transform Resumen El propósito de esta tesis es desarrollar esquemas eficientes para clasificar imágenes n-dimen- sionales usando técnicas de segmentación y morfología matemática (MM), y máquinas de soporte vectorial (SVM) como clasificadores. Por esquemas eficientes entendemos aquellos que producen buenos resultados en términos de precisión así como los que se pueden ejecutar en tiempo real en arquitecturas de bajo coste. En la búsqueda de estos esquemas establecemos pues un doble objetivo; uno, en referencia a la precisión y, otro, al tiempo de ejecución. El primer de ellos se logra mediante el diseño de nuevos esquemas, mientras que, el segundo, se consigue mediante el desarrollo de técnicas que permiten su ejecución eficiente en hardware disponible en ordenadores personales, como las CPUs multi-hilo y las unidades de procesador gráfico (GPU, Graphics Processing Unit, en inglés). Las imágenes n-dimensionales engloban tanto a las imágenes de dos y tres dimensiones, ejemplo de ello son las utilizadas en el ámbito médico, como también a las imágenes desde diez hasta cientos de dimensiones, como las imágenes multi- e hiperespectrales adquiridas en teledetección. En el análisis de imágenes multi- e hiperespectrales, un pixel se representa como un vector de valores espectrales o características al que llamamos pixel vector. Las imágenes hiperespectrales se adquieren mediante sensores ópticos que capturan diferentes longitudes de onda a la vez, desde el espectro visible hasta cerca del infrarrojo. Esta colección de datos se puede ver como un cubo hiperespectral formado por varias bandas espectrales. [66]. El primer sensor multiespectral a bordo de un satélite, el Landsat-1 en 1972, era capaz de recoger 4 bandas espectrales con una resolución espacial de 80 metros por píxel [95]. Los cientos de bandas hiperespectrales llegaron en noviembre del año 2000 con el espectrómetro Hyperion [123]. Desde entonces, la teledetección ha sido un área de investigación muy activa para el mapeo de minerales [92], la identificación de zonas urbanas; por ejemplo, para la detección de cambios urbanísticos [107], y para el análisis en la degradación de los bosques, xx Resumen entre otras [66, 128]. En este trabajo se han utilizado imágenes hiperespectrales capturadas por el Airborne Visible-infrared Imaging Spectrometer (AVIRIS) [71] y el Reflective Optics System Imaging Spectrometer (ROSIS-03) [113]. Estos dos sensores hiperespectrales cubren una gama de 0,4a0,86 µm (ROSIS-03) y de 0,4a2,4µm (AVIRIS) utilizando 115 y 224 canales espectrales, respectivamente, con una resolución espacial que varía entre 1 y 20 m / pixel. Las características especiales de las imágenes hiperespectrales, que proporcionan información detallada para cada píxel, permiten distinguir entre materiales físicos y objetos incluso a nivel de un píxel. La alta dimensionalidad de los datos presenta nuevos desafíos en técnicas como desmezclado de información espectral (unmixing) [89, 23, 24], detección de objetivos y anomalías [107, 14], extracción de características [14], o caracterización y clasificación de la superficie terrestre [96, 128]. Este último ha sido un campo muy investigado en las últimas décadas. Entre los diferentes retos, esta tesis se ocupa de la clasificación de imágenes n-dimensionales. La clasificación consiste en agrupar y etiquetar elementos que tengan en común una o más propiedades o características. El éxito de una clasificación depende de la habilidad para categorizar y/o discrimar cosas, así como en el conjunto del características usadas para tal fin; por lo tanto, usar las características más relevantes es un requisito imprescindible. En el análisis de imágenes, la clasificación es un método usado con regularidad para extraer información en medicina, vigilancia, manufacturación y en teledetección [128, 146]. Las características que intervienen en el proceso de aprendizaje para clasificar un píxel en una imagen están principalmente relacionadas con la intensidad y el color de dicho pixel, por ejemplo, los canales rojo, verde y azul en una imagen en RGB. En las imágenes hiperespectrales tenemos muchas más características que podemos usar en este proceso; sin embargo, dicho conjunto de datos espectrales resulta en ocasiones redundante, por lo que se usan técnicas para extracción de características como el análisis en componentes principales (PCA) con el objetivo de reducir la dimensionalidad y extraer las principales características que representan a estas imágenes [146]. Hemos investigado diferentes técnicas para extraer características como análisis de componentes independientes (ICA), fracción mínima de ruido (Minimum Noise Fraction en inglés), DAFE (Discriminant Analysis Feature Extraction), DBFE (Decision Boundary Feature Extraction) y NWFE (Nonparametric Weighted Feature Extraction). xxi Para sacar el máximo provecho a las imágenes n-dimensionales, se han propuesto un gran número de métodos de clasificación [72, 55, 146], como máxima verosimilitud, redes neuronales, árboles de decisión y métodos basados en máquinas de soporte vectorial (Support Vector Machines (SVM) en inglés). En un estudio exhaustivo de clasificadores presentado en [60] los autores concluyeron que las SVMs estaban entre los mejores métodos. En el campo de teledetección, SVM ha demostrado que puede obtener buenos resultados en clasificación de imágenes hiperespectrales, incluso cuando el número de muestras de entrenamiento es pequeño [72, 55]; por estos motivos, hemos usado SVM como método de clasificación en este trabajo. Independientemente de la robustez de las SVMs, estos métodos procesan cada pixel de la imagen de forma independiente considerando únicamente la información espectral; sin embargo, se ha verificado que la información espacial que se puede extraer de una imagen hiperespectral mejora la precisión de la clasificación cuando se incorpora dentro de un esquema espectral-espacial [164, 15, 54, 154, 155, 41, 19, 29, 56, 57]. Por lo tanto, los esquemas de clasificación de imágenes hiperespectrales han cambiado de clasificadores a nivel de pixel hacia esquemas de clasificación espectral-espacial. Un estudio de estos avances se recoge en [58, 24, 65]. En particular, en esta tesis estamos interesados en los esquemas que extraen la información espacial en base a técnicas de segmentación y perfiles morfológicos. Las técnicas de segmentación han sido investigadas en teledetección para extraer información de las imágenes hiperespectrales e incorporarla en los esquemas de clasificación [117, 12, 154, 155, 157, 132, 160]. Al usar las regiones creadas mediante segmentación tenemos en cuenta las estructuras espaciales que pueden estar presentes en la imagen, por ejemplo, las técnicas de segmentación usadas en [12] están basadas en métodos de conjunto de nivel (level-set en inglés). Algoritmos basados en clustering (creación de particiones en la imagen agrupando píxeles similares) se han usado en [154], y técnicas basadas en evolución de autómatas celulares fueron diseñadas y aplicadas en [132] para segmentación de imágenes hiperespectrales sin reducir la dimensionalidad de los datos. La transformada watershed se ha aplicado en los esquemas espectrales-espaciales presentados en [117, 155, 156, 158, 78]. Tilton [160] propuso un algoritmo de segmentación jerárquica denomidado HSEG (hierarchical segmentation) que combina dos métodos de segmentación para unir regiones y mejorar los resultados. Entre todas las técnicas, la transformada watershed [168] ha despertado mayor interés en los esquemas de clasificación basados en la segmentación, a pesar de que no se puede aplicar directamente a una imagen n-dimensional. El enfoque más común para la aplicación xxii Resumen de esta técnica en imágenes hiperespectrales consiste en reducir el número de dimensiones espectrales a uno, por ejemplo, mediante reducción de características (PCA) o algoritmos de gradiente vectorial (RCMG) [155]. En esta tesis hemos desarrollado un esquema de clasificación basado en segmentación (CA–WSHED–MV) [137] usando autómatas celulares para calcular de forma eficiente la transformada watershed (CA–Watershed). Este algoritmo, CA– Watershed [139, 135, 138], es apropiado para arquitecturas multi-hilo como las CPUs y las unidades de procesador gráfico (GPUs), ya que el autómata se puede dividir en bloques que se actualizan de forma independiente. La morfología matemática (MM) se define como la teoría para el análisis de estructuras espaciales [152] y se ha usado con éxito en clasificación de imágenes n-dimensionales en el campo de la teledetección, a través de perfiles morfológicos (Morphological Profiles (MP) en inglés) [124] y perfiles morfológicos basados en atributos (attribute profiles) [41] con los que es posible analizar diferentes tipos de estructuras. Los perfiles morfológicos eliminan objetos que no encajan dentro de la estructura de un elemento de análisis, llamado structuring element (SE). Aplicando sucesivamente operaciones morfológicas con elementos estructurales de distinto tamaño se crea una representación de la imagen con varios niveles de detalle. El uso en imágenes n-dimensionales es posible aplicando los perfiles a cada banda (o a un conjunto representativo) a través de perfiles morfológicos extendidos (EMPs) [121, 15], y de perfiles basados en atributos extendidos (EAPs) [41]. Los EAPs se pueden crear con diferentes atributos, por ejemplo, el área o la desviación estandard del color de una región, extrayendo más información espacial y dando lugar a los perfiles morfológicos basados en multi-atributos (EMAPs) [42]. En general, los esquemas de clasificación basados en MM usando perfiles extendidos (EMP, EAP, EMAP) han demostrado ser más eficientes en términos de precisión que los esquemas basados en segmentación [58, 65], a cambio de incrementar el coste computacional. En esta tesis hemos desarrollado un esquema de clasificación espectralespacial basado en EMP (WT–EMP), teniendo en cuenta su posterior ejecución en hardware de bajo coste para procesamiento en tiempo real. Este esquema crea un EMP a partir de un conjunto representativo de los datos y los concatena con la información espectral creando un nuevo conjunto nuevo de características para cada pixel. Independientemente de las técnicas utilizadas en los esquemas de clasificación, fallos en la calibración de los sensores o fenómenos atmosféricos pueden afectar a la calidad de los datos [146]. La consecuencia más común es la presencia de ruido en la imagen. Por lo tanto, normalmente se requiere un preprocesado de los datos para reducción de ruido [164, 166] xxiii o corrección de dispersión [140]. La transformada wavelet es una herramienta matemática para procesamiento de señales que se ha aplicado en teledetección para filtrado de ruido, así como otros preprocesados de la imagen como compresión de datos [61] y reducción de características [86]. En el esquema WT–EMP propuesto en esta tesis, usamos la transformada wavelet para extracción de características antes de crear el perfil morfológico extendido, y también para eliminar ruido en cada banda espectral de la imagen original. A pesar de que existen una gran variedad de esquemas de clasificación espectral-espacial, la mayoría resultan computacionalmente poco eficientes en términos de tiempos de ejecución debido al gran volumen de datos al que deben hacer frente. Por lo tanto, para el uso de estos esquemas en aplicaciones en tiempo real necesitamos implementaciones eficientes en las arquitecturas usadas para su ejecución. Esto es de especial interés en aplicaciones con tiempo de respuesta crítico como monitorización de desastres naturales o detección de objetivos a bordo para salvamento marítimo, donde la toma de decisiones se hace en tiempo real [129, 77, 18]. Además, como las imágenes n-dimensionales han estado más disponibles en los últimos años debido principalmente a la reducción en tamaño y coste de los sensores espectrales, el número de aplicaciones a mediana y pequeña escala como es el control de la calidad en los alimentos [30], la detección de falsificaciones en obras de arte [100], el diagnóstico de enfermedades [30, 102] y medicina forense [50], también se ha visto incrementado. La necesidad de una computación eficiente en arquitecturas de bajo coste está aumentando a medida que incorporamos esta tecnología en aplicaciones de uso diario. En esta tesis nos hemos centrado en diseñar y desarrollar esquemas de clasificación eficientes aplicados al campo de la teledetección para procesar imágenes de la superficie terrestre en tiempo real. El campo de la computación de altas prestaciones aplicado a la teledetección abarca desde grandes infraestructuras de servidores [49, 75, 130] hasta hardware de bajo coste como las FPGAs (Field Programmable Gate Arrays, en inglés) [129, 126, 67], pasando por CPUs multihilo y las GPUs [126, 127, 67, 36, 20, 137]. Mucha de la investigación en este área se ha centrado en el campo de desmezclado de información espectral [129, 67] y en la detección de objetivos [77, 18]. El hardware más idóneo depende principalmente de la aplicación final, el presupuesto y el espacio disponible para instalar dicho hardware. En el caso de esquemas de clasificación espectral-espacial muy pocos se han adaptado para la GPU [20, 137] o han sido especialmente diseñados desde el principio para tal fin. Dadas la complejidad de los sistemas de clasificación espectral-espacial y la gran cantidad de datos disponibles en las imágenes hiperespectrales, nos planteamos la siguiente pregunta xxiv Resumen en esta tesis: ¿Es posible diseñar esquemas de clasificación espectral-espacial que produzcan buenos resultados en términos de precisión, y que se puedan ejecutar en tiempo real utilizando infraestructuras de computación de bajo coste para el procesamiento a bordo de datos hiperespectrales?. El éxito de un esquema eficiente requerirá de un estudio a fondo para encontrar técnicas que mejoren la precisión de la clasificación adaptadas al modelo de computación de este hardware de bajo coste, como una GPU. Las contribuciones principales de esa tesis son: 1. Análisis de esquemas de clasificación espectral-espacial basados en la segmentación y perfiles morfológicos. En particular, nos hemos centrado en distintas formas de incorporar la información espacial en los sistemas de clasificación hiperespectrales basados en SVM. Dedicamos especial atención a técnicas para la extracción de características, así como al procesamiento espacial basado en morfología matemática, técnicas de segmentación como la transformada watershed, y técnicas de fusión de datos para combinar la información espectral y espacial. 2. Propuesta de esquemas de clasificación espectral-espacial. Hemos propuesto los siguientes esquemas: – CA-WSHED–MV [137, 136] es un esquema basado en segmentación, SVM y fusión de datos vía votación mayoritaria. Este esquema está basado en el framework de clasificación espectral-espacial propuesto por Tarabalka [156, 155]. El esquema consiste en calcular un gradiente vectorial (RCMG) que reduce la dimensionalidad de la imagen a una sola banda, la cual es segmentada usando una transformada watershed basada en autómatas celulares. La clasificación se lleva a cabo mediante SVM. Por último, las regiones segmentadas se combinan con los resultados de clasificación usando una técnica de votación. La novedad de este esquema es que el algoritmo de watershed, basado en autómatas celulares, no genera líneas de segmentación las cuales no están asignadas a ninguna región y, por lo tanto, no necesita un procesamiento adicional para incluir dichas líneas en el esquema. Además, este algoritmo sigue un modelo de computación en el cual las xxv celdas del autómata se pueden dividir en grupos y asignarse a diferentes unidades de computación, donde se pueden actualizar de forma asíncrona. – WT–EMP [133] es un esquema basado en wavelets, morfología matemática y SVM. Fue diseñado teniendo en cuenta su posterior ejecución en GPU. La transformada wavelet se usa 1) para extraer información espectral usando filtros 9/7, y 2) para reducción de ruido utilizando un conjunto de tres filtros [148]. La morfología matemática se usa para crear un perfil morfológico extendido a partir de las características extraídas por wavelets. Este nuevo esquema de clasificación espectral-espacial mejora los resultados en términos de clasificación en comparación con otros basados en segmentación y MM. 3. Desarrollo de técnicas y estrategias para computación eficiente en GPU. Hemos aplicado diferentes estrategias, y en particular, hemos propuesto una estrategia basada en computación asíncrona (block–asynchronous) con la que es posible ejecutar autómatas celulares en GPU de forma más eficiente. Esta estrategia reduce el número de puntos globales de sincronización y explota de forma eficiente la jerarquía de memoria de esta arquitectura; además, resulta adecuada para arquitecturas multi-hilo y se ha adaptado para su uso tanto en imágenes 2D como en volúmenes en 3D. En particular, la hemos aplicado a: – CA–Watershed: éste es el autómata celular asíncrono para calcular la transformada watershed en GPU [139, 135, 138]. El comportamiento asíncrono de esta propuesta introduce irregularidades en los bordes entre regiones que hemos solucionado corrigiendo la velocidad de propagación entre bloques, usando una técnica denominada wavefront [112]. – Operaciones morfológicas de apertura y cierre por reconstrucción basados en el mismo principio de computación asíncrona en GPU [134]. Estas técnicas se usan para crear los perfiles morfológicos empleados en el esquema WT-EMP. – Filtrado basado en atributos. Se trata de una nueva propuesta para el filtrado (apertura y cierre) de imágenes en escala de grises en GPU y que además, se puede extender a diferentes atributos. Este tipo de filtrado se aplica en esquemas de clasificación basados en atributos extendidos (EAPs) y multi-atributos (EMAPs). 4. Implementación eficiente de esquemas de clasificación espectral-espacial en CPUs multi-hilo y GPUs, usando OpenMP y CUDa, respectivamente. 4Chapter 1. Thesis overview of Morphological Profiles (MPs) [124] and Attribute Profiles (APs) [41], which are used to model different kinds of spatial structures (objects). The profiles keep objects in the image if a structuring element or attribute fits within the object, otherwise they are removed. Through the sequential application of morphological operations and increasing the size of the structuring element (or filter), the MP and AP create a multilevel characterization of the image. The extension to n-dimensional images is possible by applying the profile to each band (or to a representative subset) through the Extended Morphological Profile (EMP) [121, 15] and the Extended Attribute Profile (EAP) [41], respectively. The latter can be further extended with different attribute filters, (e.g. the area and standard deviation of pixel colors within a region), extracting more spatial information from the image, creating an Extended Multi-Attribute Profile (EMAP) [42]. In general, MM-based classification schemes using extended profiles have shown to be more efficient than segmentation-based approaches in terms of classification accuracy [58, 65], at the expense of increasing the computational cost. In this thesis we have developed a scheme based on Extended Morphological Profiles (WT–EMP) and taking into account its subsequent projection on low-cost computing infrastructures for real-time processing. The scheme creates the EMP from a representative subset of the hyperspectral bands that is combined with the spectral data in a new vector of features. Regardless of the techniques used in remote sensing, undesirable artifacts such as improper calibration of the sensor or atmospheric phenomena may affect data quality [146]. The most common consequence of these artifacts is the presence of noise in the image. Therefore, a preprocessing step, such as denoising [164, 166] or scatter correction [140], is usually required before classifying the images. Wavelets are mathematical tools for signal processing and have been investigated in remote sensing for filtering the noise introduced in the acquisition of the image, as well as additional preprocessing like data compression [61] and feature extraction [86]. In the WT–EMP scheme proposed in this thesis, wavelets are used for feature extraction to create the EMP, and also for denoising over each band of the original hyperspectral image. Although there are many spectral-spatial classification schemes, most are computationally inefficient in terms of execution time, as they have to deal with a high number of features resulting in large execution times. Therefore, the computation of these schemes for real-time applications requires their efficient implementation on the adequate computing architecture. This is particularly true in time-critical applications in remote sensing, such as natural disasters monitoring and on-board target detection for maritime rescue where decisions are made 1.1. Main contributions 5 in real time [129, 77, 18]. In addition, as the hyperspectral images have been widely available in recent years owing to the reduction in the size and cost of the sensors, the number of applications at lab scale, such as food quality control [30], art forgery detection [100], disease diagnosis [30, 102] and forensics [50], has also increased. The need for efficient computation on low-cost computing infrastructures is increasing in line with the incorporation of technology into everyday applications. This thesis focuses on developing efficient spectral-spatial schemes for land-cover classification in remote sensing imaging for real-time applications. The research in High Performance Computing (HPC) for remote sensing applications covers a field ranging from infrastructures of clusters [49, 75, 130] to Field Programmable Gate Arrays (FPGAs) [129, 126, 67] and commodity hardware such as multi-threaded CPUs and many-core GPUs [126, 127, 67, 36, 20, 137]. Most of the research in HPC for remote sensing has been done in the field of spectral unmixing involving endmember extraction [129, 67], and also in target detection [77, 18]. The most suitable hardware for HPC depends mainly on the final application, the budget and the space available for the computing infrastructures. The GPU, with its high computational capacity has not been yet fully exploited. In the case of spectral-spatial classification schemes, only a few have been adapted for GPU [20, 137], and even fewer have been specially designed for that purpose from the outset. Given the complexity of the spectral-spatial classification schemes, and the high amount of data available in the hyperspectral images, we posed the following question in this thesis: Is it possible to design efficient spectral-spatial classification schemes that produce good classification results and can be executed in real-time, using low-cost computing infrastructures for on-board processing of hyperspectral information? The success of an efficient scheme will require a thorough study to find techniques that improve the classification accuracy, while matching the computing model of the low-cost computing infrastructures, such as a GPU. 1.1 Main contributions As a result of the research conducted in this thesis to find a solution to our main question, the following contributions to the field of remote sensing and HPC have been produced: 6Chapter 1. Thesis overview 1. Analysis of spectral-spatial classification schemes based on segmentation and morphological profiles. In particular, we focus on the different ways of incorporating spatial information into the pixel-wise spectral classification schemes based on the SVM classifier. We devote special attention to feature extraction techniques, as well as spatial processing by mathematical morphology, segmentation techniques based on clustering, such as kmeans and quick-shift, and segmentation techniques based on region growing such as the watershed transform, and data fusion strategies for combining the spectral and spatial information. We have taken into account the efficient computation on GPU of the techniques under study. 2. Proposal of spectral-spatial classification schemes. The following schemes are proposed: – CA-WSHED–MV [137, 136] is a scheme based on segmentation, SVM and Majority Vote. This scheme is based on the spectral-spatial classification framework proposed by Tarabalka et al. [156, 155]. The scheme consists of the calculation of a Robust Color Morphological Gradient (RCMG) which reduces the dimensionality of the hyperspectral image, followed by the calculation of a watershed transform based on cellular automata which produces the spatial results. The classification is carried out by SVM. Finally, the spectral and spatial results are combined with a majority vote. The novelty of this scheme is mainly introduced in the watershed algorithm based on cellular automata used for segmenting the hyperspectral image. This algorithm follows a computational model in which the grid of cells of the automaton is partitioned into regular regions that are assigned to different blocks of threads on the GPU that can be asynchronously updated. In addition, this implementation does not create the so-called watershed lines, and thus it is not necessary to compute a standard vector median [7] for every watershed region as in [155]. – WT–EMP [133] is a scheme based on wavelets, MM and SVM. This scheme was designed taking into account its efficient computation in a subsequent GPU execution. The Wavelet Transform is used for feature extraction and image denoising using 9/7 wavelet filters for the former and a set of three filters for perfect reconstruction [148] for the latter. Mathematical Morphology is used for creating the EMP from the features extracted by wavelets. This new spectral-spatial classifi- 1.1. Main contributions 7 cation scheme improves the classification results in terms of accuracy in comparison to other spectral-spatial schemes based on segmentation and Mathematical Morphology. 3. Development of techniques and strategies for efficient GPU computing. Different strategies are applied and, in particular, a block–asynchronous strategy that maps cellular automata on the GPU is proposed. This strategy reduces the number of points of global synchronization allowing efficient exploitation of the memory hierarchy of this architecture. The block–asynchronous strategy is also adequate for multicore architectures, and it has been tuned to be used with 2D and 3D images. It is applied to: – CA–Watershed: this is an asynchronous cellular automaton to compute the watershed transform on GPU [139, 135, 138]. The asynchronous behavior of the CA–Watershed introduces artifacts in the border of the segmented regions that are fixed by correcting the data propagation speed among the blocks using wavefront techniques [112]. – Opening and closing: an improved algorithm for opening and closing by reconstruction on GPU based on the block-asynchronous approach (BAR) [134]. The algorithm is applied to create the Extended Morphological Profile used in landcover spectral-spatial classification schemes. – Attribute filtering: a new proposal for greyscale attribute opening and closing on GPU. This proposal is the first attempt to compute this filter on greyscale images on GPU and can be extended to other attributes. The proposal can be applied to create an Extended Attribute Profile for land-cover spectral-spatial classification schemes. 4. Efficient implementation of spectral-spatial classification schemes on multi-threaded CPUs and many-core GPUs by using OpenMP and CUDA, respectively. – CA–WSHED-GPU [136, 137] is an efficient projection on GPU of the first proposed scheme (CA-WSHED-MV). Different hyperspectral data partitioning strategies and thread block arrangements are studied in order to effectively exploit the memory and computing capabilities of this architecture. The spectral partitioning is used to compute the RCMG. The watershed transform is based on the asynchronous cellular automaton (CA-Watershed) that efficiently exploits the GPU 8Chapter 1. Thesis overview architecture, and the majority vote is done by atomic operations to avoid memory race conditions. – WT–EMP–GPU [134] is a GPU implementation of the WT–EMP scheme that achieves real-time in commodity hardware. We adapted the feature extraction by wavelets computing thousands of pixel vectors in parallel. A new implementation for two dimensional wavelet transforms was required in order to manage the three filters used in the denoising step. Finally, we used the block–asynchronous strategy for morphological reconstruction to compute the EMP used in the scheme. Despite the different GPU solutions for SVM classification found in the literature, a new implementation was developed for classifying multi-class problems on GPU [134], that is compatible with the one against one trained models produced by the facto LIBSVM library [35]. 1.2 Publications 1.2.1 Book Chapters [137] P. Quesada-Barriuso, F. Argüello, and D. B. Heras, “Computing efficiently spectralspatial classification of hyperspectral images on commodity gpus,” in Recent Advances in Knowledge-based Paradigms and Applications (J. W. Tweedale and L. C. Jain, eds.), vol. 234 of Advances in Intelligent Systems and Computing, Ch. 2, pp. 19–42, Springer International Publishing, 2014. 1.2.2 International Journals [138] P. Quesada-Barriuso, D. B. Heras, and F. Argüello, “Efficient 2D and 3D watershed on graphics processing unit: block-asynchronous approaches based on cellular automata,” Computers & Electrical Engineering, vol. 39, no. 8, pp. 2638–2655, 2013. [78] D. B. Heras, F. Argüello, and P. Quesada-Barriuso, “Exploring ELM-based spatialspectral classification of hyperspectral images,” International Journal of Remote Sensing, vol. 35, no. 2, pp. 401–423, 2014. [133] P. Quesada-Barriuso, F. Argüello, and D. B. Heras, “Spectral-spatial classification of hyperspectral images using wavelets and extended morphological profiles,” Selected Topics 1.2. Publications 9 in Applied Earth Observations and Remote Sensing, IEEE Journal of, vol. 7, no. 4, pp. 1177–1185, 2014. [134] P. Quesada-Barriuso, F. Argüello, D. B. Heras, and J. A. Benediktsson, “Wavelet-based classification of hyperspectral images using extended morphological profiles on graphics processing units,” Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, vol. PP, no. 99, pp. 1–9, 2015 (published online, print edition pending). [101] J. López-Fandiño, P. Quesada-Barriuso, D. B. Heras, and F. Argüello, “Efficient ELM- based techniques for the classification of hyperspectral remote sensing images on commodity gpus,” Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, vol. PP, no. 99, pp. 1–10, 2015 (published online, print edition pending). 1.2.3 International Conferences [135] P. Quesada-Barriuso, D. B. Heras, and F. Argüello, “Efficient GPU asynchronous implementation of a watershed algorithm based on cellular automata,” in Parallel and Distributed Processing with Applications (ISPA), 2012 IEEE 10th International Symposium on, pp. 79– 86, 2012. [136] P. Quesada-Barriuso, F. Argüello, and D. B. Heras, “Efficient segmentation of hyperspectral images on commodity GPUs,” in 16th International Conference on Knowledge-Based and Intelligent Information & Engineering System, vol. 243, pp. 2130–2139, 2012. 1.2.4 National Conferences [139] P. Quesada-Barriuso, J. Lamas-Rodríguez, D. B. Heras, and F. Argüello, “Influencia de las mesetas en la implementación de watershed sobre GPUs,” in XXIII Jornadas de Paralelismo, pp. 249–254, 2012. [94] J. Lamas-Rodríguez, P. Quesada-Barriuso, F. Argüello, D. B. Heras and y M. Bóo, “Proyección del método de segmentación del conjunto de nivel en GPU,” in XXIII Jornadas de Paralelismo, pp. 273–278, 2012. 10 Chapter 1. Thesis overview 1.3 Thesis organization This thesis is organized in six chapters. Chapter 1 introduces the objectives and the main contributions of this thesis. Chapter 2 presents the fundamental concepts involved in this work. First, the main characteristics of the n-dimensional images are reviewed. Second, the techniques used in the proposed spectral-spatial classification schemes are detailed. Different feature extraction techniques that have been successfully applied in remote sensing classification are described. The concepts regarding the spatial processing of the images, such as the segmentation, mathematical morphology and attribute filtering are also reviewed in this chapter. Then, cellular automata and the wavelet transform are introduced, and we review the pixel-wise classifiers, in particular the SVM, the one used in our schemes. In addition, the OpenMP and Compute Unified Device Architecture (CUDA) computing model, as well as the hardware resources and the performance measures used in the experiments are introduced in this chapter. Finally, the datasets used in the experiments are described. Chapter 3 presents our two schemes proposed for efficient n-dimensional image classification: CA–WSHED–MV and WT–EMP. The general framework of the spectral-spatial classification schemes and the common data fusion strategies to join the spectral and the spatial information of such schemes are presented in this chapter. The results thereof are compared on real hyperspectral images, in terms of classification accuracy, with the classification schemes found in the literature based on segmentation and MM. Chapter 4 focuses on the techniques and strategies applied for the efficient computation of the proposed schemes on commodity hardware such as multi-threaded CPUs and many-core GPUs. First, we describe the general strategies related to data partitioning, data movement and data packing, and highlight the challenges of GPU computing. Then, we present the block–asynchronous computation strategy proposed in this thesis. This approach is applied to compute an asynchronous CA–Watershed, opening and closing by reconstruction and attribute filtering by opening and closing. The strategies for efficient projection of wavelets, as well as multi-class SVM classification on GPU are also described in this chapter. The techniques are independently analyzed comparing the execution times to parallel multi-threaded CPU implementations using OpenMP. In chapter 5 we describe the CUDA implementation of the proposed schemes, named CA–WSHED–GPU and WT–EMP–GPU, applying the techniques described in chapter 4. We also analyze the optimal GPU parameters and the hardware resources that usually limit the occupancy on this architecture. The performance is measured in terms of execution time 1.3. Thesis organization 11 and speedup on real hyperspectral images over CPU multi-threaded implementations using commodity hardware. Finally, the last chapter remarks the contributions, presents the conclusions and proposes future research developments. CHAPTER 2 FUNDAMENTALS In this chapter we review the fundamental concepts and present the techniques used in the spectral-spatial classification schemes proposed in this thesis. In particular, the main characteristics of the n-dimensional images are reviewed and the concepts regarding the spatial processing of the images, such as the watershed transform, cellular automata, vectorial gradients, mathematical morphology and wavelet transform are outlined. We then go on to describe the pixel-wise classifiers and the classification map produced by the same, in particular for SVM that is the classifier used in our scheme. In addition, the OpenMP and the CUDA parallel programming models, as well as the hardware resources and the performance measures are introduced in this chapter. Finally, the datasets used in the experiments are detailed. 2.1 n-dimensional images An-dimensional (nD) image is considered a dataset where two dimensions represent a spatial location [x,y], while the remaining dimensions represent phenomena occurring per spatial location [3]. Therefore, a pixel is represented as a vector of values or features (pixel vector). For example, a nD image can be produced by optical sensors which capture different wavelength responses (one per spectral band) in the same line scan. Multispectral images usually are composed of a few, usually three to ten, spectral channels, while hyperspectral images have hundreds of spectral bands. These datasets can be seen as spectral cubes of several (hundreds) images where [x,y]varies across the physical space and [z]varies across the electromagnetic spectrum. [66]. A nD image can also be a set of images of the same scene acquired by the 20 Chapter 2. Fundamentals Input: input data X Output: partitioning of the data in kclusters 1: choose the number of kclusters and initialize their position cjto kfeature vectors 2: repeat 3: for each feature vector x∈X do 4: assign xto its closest cluster based on the Euclidean distance 5: end for 6: for each cluster Cjdo 7: recompute cjto the mean of all the points x(j) ibelonging to that cluster 8: end for 9: compute the square error function (2.10) 10: until stability Figure 2.2: Pseudocode for k-means clustering. Input: input image X Output: hierarchical image segmentation 1: label each pixel of Xas a separate region 2: while the number of regions remaining is greater than two do 3: repeat 4: calculate the dissimilarity criterion di,jbetween each spatially adjacent region 5: find the min(di,j)and merge all pairs of spatially adjacent regions with meet that value 6: calculate the dissimilarity criterion si,jbetween all pairs of non-spatially adjacent regions 7: merge all pairs of non-spatially adjacent regions with si,j≤min(di,j) 8: until stability 9: save current segmentation map 10: end while Figure 2.3: Pseudocode for HSEG algorithm. By analogy, hierarchical segmentation can be defined as a family of fine to coarse image partitions [153]. Tilton proposed a Hierarchical Image Segmentation (HSEG) algorithm [161] for exploiting the information content on the segmentation hierarchy. It is a hybrid segmentation technique based on hierarchical step-wise optimization (HSWO) [13] and data clustering. The algorithm alternately performs region growing merging spatially adjacent regions, and spectral clustering merging spatially non-adjacent regions. The HSEG pseudocode shown in Figure 2.3 is based on the description given in [161, 173]. The inner loop (lines 3–8) merges the segmented regions of the image. Lines 4 and 2.3. Segmentation techniques 21 5 correspond to the region growing step and lines 6 and 7 constitute the spectral clustering part. The dissimilarity criterion is based on the Euclidean distance characterized by the mean vectors of the regions. Stability is reached when the number of regions remaining in the segmentation map is lower than a preset value of regions [161]. The outer loop, lines 2–10 performs the hierarchical segmentation of the image. In clustering-based image segmentation, the labels in the clustering map are created from kdifferent classes; i.e., there are only kdistinct values and the segmented regions are not connected by a unique label. Therefore, a Connected Component Labelling (CCL) algorithm is required to re-label the clusters and produce a segmentation map [83]. The segmentation of nD images aims at finding distinct structures in the spectral domain. 2.3.2 Watershed transform The watershed transform is a widely used method for non-supervised image segmentation, especially suitable for low-contrast images. The idea behind this method comes from geography. A greyscale image can be represented as a topographic relief, where the height of each pixel is directly related to its grey level. The dividing lines of the catchment basins for precipitation falling over the region are called watershed lines [168]. Various definitions, algorithms and implementations can be found in the literature but, in practice, they can be classified into two groups: those based on the specification of a recursive algorithm by Vincent and Soille [168], and those based on the distance functions defined by Meyer [111]. An intuitive approach is to imagine the terrain being immersed in a lake, with holes pierced in local minima [142]. By flooding the terrain into the water, catchment basins will fill up with water as illustrated in Figure 2.4. At the points where water coming from different basins would meet, dams are built. The process stops when the water level reaches the highest peak. The terrain is partitioned into regions separated by the watershed lines. In Figure 2.4 the image has two minimum grey values representing two valleys in the terrain. One of the main advantages of the watershed transform is that all regions of the image are well defined at the end of the segmentation process, even if the contrast of the image is poor. Hence it has been widely used in image processing, biomedicine and physics. Nonetheless, the results are over-segmented owing to the large number of regions detected. This problem is overcome by preprocessing the image with the objective of reducing the number of regions, for example by a marked-controlled watershed transform that preselects the regions of interest [22]. 22 Chapter 2. Fundamentals Figure 2.4: Flooding process in a one dimensional image with two minimum, generating a watershed line at the end of the process. Let us introduce a few concepts and notations about topography in order to continue with the watershed transform. A greyscale image may be considered as a graph G= (V,A)with a finite set of Vvertexes (pixels) and a set of arcs A⊆V×Vdefining the connectivity. Two pixels uand vare connected if (u,v)∈A. The pixels connected to u, called neighbors, are denoted by N(u). The most widely used connectivity is four, considering the orthogonal neighbors, left, right, up and down, known as Von Neumann connectivity. Another variation is the Moore neighborhood, where the eight neighbors surrounding a pixel are connected. Figure 2.5(a) shows a pixel with 4- and 8-connectivity. The slope between two neighbors is defined by: ∀u∈V,∀v∈N(u),slope(u,v) = h(u)−h(v), where h(u)is the grey value (altitude) of the pixel u. The lower slope is defined as the maximal slope connecting uto any of its neighbors of a lower altitude: LS(u) = max(h(u)−h(v))|v∈Γ(u), with Γ(u)the set of neighbors vwith h(v)<h(u). If Γ(u)has more than one element, NF(u) represents an arbitrary element of that set. (a) (b) Figure 2.5: (a) Pixel 4– and 8–connectivity, (b) backward N+ G(p)and forward N− G(p)neighborhood of a pixel p with 8–connectivity. 2.3. Segmentation techniques 23 (a) (b) (c) (d) Figure 2.6: Watershed based on Hill-Climbing algorithm. Example of 1D image represented as a terrain by dashed lines and as grey values by the squares: (a) detecting and labelling all minima in the image with unique labels (“A” and “B”), (b) and (c) continues propagating the labels upwards, climbing up the hill, (d) result of the segmentation in two regions. A plateau is a connected subgraph P= (VP,AP)⊆G, where ∀u∈VP,h(u) = c, and cis the altitude of the plateau. Thus, a plateau is a region of constant grey value within the image. The set N=(u)denotes the pixels v∈N(u)with h(u) = h(v). A connected component of a graph Gis a subgraph Pin which pixels are connected to each other by paths. A plateau Pis a level component of the image, considered as a valued graph, i.e., a connected component of pixels of constant grey value h[142]. Finally, the lower border of a plateau Pis defined as: ∂− P={u∈VP|∃v∈N(u),h(v)<h(u)}. If ∂− P=/0 the plateau is called a minimum plateau. All the pixels within a minimum plateau are also minimum. In contrast, if ∂− P6=/0 the plateau is called a non-minimum plateau. Different algorithms have been implemented to compute the watershed transform using sequential structures as queues or graphs to simulate the flooding process [142]. Although various implementations can be found in the literature, in this thesis we follow the Hill-Climbing algorithm based on the topographical distance by Meyer [111]. This algorithm starts by detecting and labelling all minima in the image with unique labels, as illustrated in Figure 2.6(a). The process continues by propagating the labels upwards, climbing up the hill, following the path defined by the lower slope of each pixel, Figure 2.6(b) and Figure 2.6(c). The result of the segmentation, as shown in Figure 2.6(d), is a set of regions, each one represented by a catchment basin, with their own label. At the end all the pixels belong to a region and the watershed lines are the limits between these regions. 24 Chapter 2. Fundamentals Problems arise for digital images with plateaus, as it is not possible to know a priori whether a plateau is minimum or non-minimum, so an additional processing is required. The most common solution is to preprocess the image by calculating its lower complete image [111], where each pixel has at least one neighbor with a lower value, except those pixels which are minima. Another alternative is to calculate the distances of an inner pixel to the lower border of the plateau during the watershed processing, which is the strategy selected in this work. The pseudocode of the watershed transform based on the Hill-Climbing algorithm is given in Section 3.2.2 where its implementation based on cellular automata is explained in detail. 2.4 Cellular automata Cellular Automata (CA) constitute a computing model that has been extensively used for artificial life [51], pattern recognition [37] or image processing [143]. The popularity of CA is mainly due to the simplicity of modelling complex problems with the help of local information only. CA are composed of a set of cells arranged into a regular grid of one, two or three dimensions originally proposed by John Von Neumann [116] as formal models of self-reproducing organisms. In the case of two dimensions, each cell is connected to its four or eight adjacent neighbors, depending on the connectivity, as shown in Fig. 2.7. The most widely used connectivity is the Von Neumann neighborhood, considering the orthogonal neighbors, left, right, up and down. The Moore neighborhood, where the eight neighbors surrounding a cell are connected is a variation also used for connecting cells of the automaton. Figure 2.7: Cells arranged in a regular grid of two dimensions. Cells are connected to four (left) and eight (right) adjacent neighbors. 2.4. Cellular automata 25 Figure 2.8: Synchronous updating process in a cellular automaton. The updates of the cells are performed synchronously and in discrete time steps. Grey squares represent cells whose state has changed. Cells in a cellular automaton can be in one of a finite number of possible states. Each cell changes its state depending on the current state and the state of its neighbors. This requires a strict order for updating the automaton, where a cell cannot be updated until all other cells have also been updated. The updates of the cells are usually performed synchronously and in discrete time steps. Figure 2.8 shows an example of a 6×6 automaton, using 8-connectivity, which is updated synchronously. As illustrated in the figure, the updating process takes places in discrete time steps. The cells that are updated in this example are shaded in grey. If the updates of the cells are not required to take place synchronously, but each one can be updated to its next state an unbounded number of times without synchronization, then we have an asynchronous automaton [115]. In this case, the grid can be partitioned into different regions which can be updated independently. It is possible to ignore the synchronization points associated to the evolution of the automaton, resulting in a so-called asynchronous automaton. In this case, however, the correctness and convergence of the algorithm could be severely affected [4]. In order to efficiently develop asynchronous computing schemes, it is important to investigate the non-deterministic and probabilistic behavior associated to such schemes [2]. Two dimensional CA can be directly mapped into 2D images where each cell in the automaton corresponds to a pixel in the image. The connectivity used in the automaton is denoted in the image by N(u). The same can be applied for 3D cellular automata and 3D images by modifying the connectivity to take into account the third dimension. 26 Chapter 2. Fundamentals (a) (b) (c) Figure 2.9: (a) Original image, (b) erosion (shrinks brighter objects), and (c) dilation (expands brighter objects). The SE was a square of 5×5 pixels. 2.5 Mathematical morphology In this section we present an overview of Mathematical Morphology (MM) operators. In particular, we describe morphological operators based on the geodesic reconstruction, which are tools defined in the MM framework [152]. The Mathematical Morphology (MM) is a theory for analyzing and processing spatial structures from images [149]. The techniques are based on two basic operators, erosion (ε) and dilation (δ), from which a set of advanced analysis tools is constructed. Erosion and dilation transform an image Iusing a structuring element (SE), giving as output for each pixel pthe infimum (V)or supremum (W), respectively, of the intensity values of the set of pixels included by the SE when it is centered on p. Generally speaking, the erosion operator shrinks objects that are brighter than their surroundings, whereas the dilation operator expands them. Figure 2.9 shows the erosion and dilation of the image of Lena using a square SE of 5×5 pixels. 2.5.1 Opening and closing by reconstruction The dilation of an eroded image is known as opening, γ(I). Conversely, the erosion of a dilated image is known as closing, φ(I). The opening operator flattens bright objects of an image if the SE fits within the objects, and the closing operator has the opposite effect. Therefore, opening and closing are used to extract information related to the shape and size of the objects. However, opening and closing are not connected filters; i.e., these operators 2.5. Mathematical morphology 27 (a) (b) (c) Figure 2.10: A subset of 256×256 pixels of the image of Lena. (a) erosion by a SE of 5×5 pixels, (b) opening, and (c) opening by reconstruction. do not preserve edges between objects as they work on pixels rather than image structures. For instance, adjacent regions can be merged into one, and thus this biases the analysis of the spatial distribution. It is possible to construct opening and closing operators based on the geodesic reconstruction to completely preserve or remove the spatial structures of an image if the SE fits (or not) within the objects [149]. The geodesic dilation of an image I(the mask) from J(the marker) is defined as: δ(1) I(J) = δ(1)(J)∧I,(2.11) where δ(1)denotes the elementary dilation and ∧the point-wise minimum. The opening by reconstruction γ(n) r(I)of an image Iis defined as the reconstruction by dilation of Ifrom the erosion with a SE of size nof I. First, the image is transformed by an erosion εn(I)creating a marker image J. Then, reconstruction by dilation Rδ Iis an iterative process that applies geodesic dilation on the marker image until stability (δ(n) I=δ(n+1) I): Rδ I(J) = δ(n) I(J) = δ(1) Ihδ(n−1) I(J)i.(2.12) The reconstruction permits the full retrieval of all those structures that were not completely removed by the erosion. Thus, the morphological reconstruction needs several iterations before stability is attained. Figure 2.10 shows the result of opening and opening by reconstruction over a subset of 256×256 pixels of the image of Lena. It can be observed in Figure 2.10(b) that the edges in the hair or in the feathers of the hat are not preserved in the opening. How- 28 Chapter 2. Fundamentals ever, in Figure 2.10(c) the structures of the hair and feathers that were not completely removed by the erosion are successfully reconstructed by opening by reconstruction. By duality, the closing by reconstruction φ(n) r(I)of an image Iis defined as the reconstruction by erosion of Ifrom the dilation with a SE of size nof I. 2.5.2 Attribute filtering Often, a classic image analysis preprocessing problem consists of filtering out small light (or dark) particles from greyscale images without damaging the remaining structures [170]. Opening and closing by reconstruction are good operators for this task. However, these operators only extract information related to the size of the objects. Therefore, if the structures to be preserved are elongated objects, they can be completely removed. Morphological attribute filters are connected operators that process an image according to a criterion [25], such as the area or the standard deviation of the pixel intensity, removing plateaus (connected components) that do not satisfy a given criterion. Morphological attribute filters are connected filters. The following definitions are for binary images, although generalization to grey-scale images can be achieved through threshold decomposition [76]. The connected opening Γxof a set Xat a point xis defined as: Γx(X) =    Xif x∈X, /0 otherwise. (2.13) Figure 2.11 shows an example of a connected opening on a binary image with a set Xconsisting of three connected components, named C1 1,C2 1, and C3 1. The i-th component is indicated by the superscript, and the background (0) or foreground (1) is indicated by the subscript. The opening Γxpreserves only that connected component in Xwhich contains the pixel x, as illustrated in Figure 2.11(b). A trivial opening ΓTuses an increasing criterion Tto filter connected components and it is defined as follows [150]: ΓT(C) =    Cif Csatisfies criterion T, /0 otherwise. (2.14) The trivial opening preserves the connected components for which the increasing criterion T holds. For example, the following is an increasing criterion: must have an area of λpixels or more [25]. 2.5. Mathematical morphology 29 (a) (b) Figure 2.11: Example of a connected opening. (a) Binary image with three connected components named C1 1,C2 1, and C3 1, and (b) the connected opening Γxindicated by the intersection of the dot lines. The attribute opening ΓTof a set Xcombines the connected opening (2.13) and the trivial opening (2.14) and is defined as follows: ΓT(X) = [ x∈X ΓT(Γx(X)).(2.15) The attribute opening is the union of all connected components of Xwhich meet the criterion T. The attribute opening can be explained by using a tree representation of the image, as illustrated in Figure 2.12(b). The root of the tree represents the background of the image, and the leaves corresponds to the connected components in the foreground. The attribute filtering, (a) (b) (c) (d) Figure 2.12: Example of a attribute opening. (a) Binary image with three connected components named C1 1,C2 1, and C3 1, (b) tree representation of the image, (c) attribute filtering removing a node of the tree, and (d) the final result. 36 Chapter 2. Fundamentals The classification phase is performed by computing the sign of the following decision function: D(x) = ∑ i∈S αiyi(xi·x)+b,(2.24) where αiare the nonzero Lagrange multipliers, Sis the subset of training samples xi, and xis the unknown sample. The parameters (αi,S,b)are found during the training phase. For each x, the corresponding class is +1 if D(x)>0 and −1 otherwise. Although the solution may be determined optimally by the SVM when the training samples are not linearly separable, the SVM approach uses kernel tricks to map the non-linearly separable data into a higher dimensional space, in order to enhance linear separability between the two classes. Several types of kernels, such as linear, polynomial, splines or Radial Basis Function (RBF) kernels can be used in SVM classifiers. We refer the reader to [1, 72, 27] for further details on SVM and kernel tricks. For hyperspectral image classification, the Gaussian RBF has been widely used. The discriminant function (2.24) can be expressed using the RBF kernel as: D(x) = ∑ i∈S αiyiexp−γ||xi−x||2+b,(2.25) where exp(−γ||xi−x||2)is the RBF kernel and γis a positive parameter controlling the width of the kernel. 2.7.1 5-fold cross validation The SVM is mainly a nonparametric method, although some parameters need to be tuned before the optimization. In the RBF kernel case, there are two parameters: C, which controls the penalty term, and γ, which is the width of the kernel. The best parameters are usually found by k–fold cross-validation, where several parameters are tested, for example by a search in a range of values (grid search). Given a set of training samples, the k-fold cross validation divides the total number of samples into kpartitions. One partition is used to validate the accuracy of the classification, and the remaining k−1 partitions are used as training data. This validation is repeated k times, using each partition once as a test. The final classification accuracy is determined by the average with the kclassification results. In this work we have used 5–fold cross validation for tuning the SVM parameters. 2.7. Pixel-wise classification by SVM 37 2.7.2 Multi-class SVM classification The standard SVM classifier is designed for two-class problems. When the data have K>2 classes a method for solving the multiclass problem must be defined. These methods may be of the type one-against-all (OAA), one-against-one (OAO) and all-at-once, where a set of binary classifiers is combined. In the OAA approach, the Kclass problem is converted into K binary classifiers where the ith classifier separates the class ifrom the remaining classes. The OAO multiclass approach creates T=K(K−1)/2 binary classifiers on each pair of classes, where Kis the number of classes. The final classification for each point is given by the highest number of votes obtained for a class in the Tclassifiers. Hsu and Lin [80] found that the OAO approach is more suitable for practical use than the other methods, mainly because the training time is shorter. This is the approach used in all the SVM classifications carried out in this work. 2.7.3 LIBSVM: the facto library for SVM In this section we present the facto library for SVM (LIBSVM) and the principal interfaces and extensions that are based on this library. LIBSVM [35] is a library for Support Vector Machine (SVM) written in C, widely used in machine learning. The library implements the One-Against-One approach for multiclass classification. It is currently one of the most widely used SVM software and it has been extended to third-party software such as Weka3, R4, Octave and Matlab5. Owing to its popularity, it has been also implemented in other languages such as Java5, Python5and C#6, as well as in other architectures, such as Cell processors and GPUs. While efficient algorithms for training SVM are available, dealing with large datasets makes training and classification a computationally challenging problem. In [109] the training phase implemented in the LIBSVM is speeded up by parallel computation in Cell processor architectures. LIBSVM has also been accelerated with GPUs using CUDA [8]. This GPU- accelerated LIBSVM implementation is a modification of the original LIBSVM that exploits CUDA with the same functionality and interface of LIBSVM. A comprehensive list of other SVM implementations and wrappers can be found in [60]. In this thesis we have implemented 3https://weka.wikispaces.com/LibSVM/ 4http://cran.r-project.org/web/packages/e1071/ 5Matlab, Octave, Java and Python implementations are included in the LIBSVM package. 6https://github.com/ccerhan/LibSVMsharp/ 38 Chapter 2. Fundamentals a new GPU proposal (GPUSVM) for the classification stage, and the LIBSVM [35] has been used as a basis for evaluating the proposed schemes in CPU. 2.8 Parallel programming models This section describes the OpenMP and CUDA programming models. First, we will describe OpenMP, an API for writing CPU multi-threaded applications on shared memory architectures. Then, we introduce the Compute Unified Device Architecture (CUDA) developed by NVIDIA, a parallel computing platform and a programming model that leverages the high computational throughput of NVIDIA GPUs. We will highlight the major differences found in the different GPU architectures used in this work. 2.8.1 OpenMP OpenMP is the standard Application Program Interface (API) for multi-threaded parallel programming on shared memory architectures [120]. Communication and coordination between threads is expressed through read/write instructions of shared variables and other synchronization mechanisms. It comprises compiler directives, library routines and environment variables and is based on a fork-join model, as illustrated in Figure 2.18 where a master thread creates a team of threads that work together in a Single Program Multiple Data (SPMD) way (different cores execute different threads operating on different data). In shared memory architectures, OpenMP threads access the same global memory where data can be shared among them or can be private for each one. From a programming perspective, data transfer for each thread is transparent and synchronization is mostly implicit (see Figure 2.18). When a thread enters a parallel region, it becomes the master, creates a thread team and forks the execution of the code among the threads and itself. At the end of the parallel region, the threads join together and the master resumes the execution of the sequential code. Different types of worksharing constructs can be used to share the work of a parallel region among the threads [33]. The loop construct distributes the iterations of one or more nested loops into chunks, among the threads in the team. By default, there is an implicit barrier at the end of a loop construct. The way the iterations are split depends on the schedule used in the loop construct [33]. On the other hand, the single construct assigns the work on only one of the threads in the team. The remaining threads wait until the end of the single construct 2.8. Parallel programming models 39 Figure 2.18: OpenMP fork-join parallel model. A master thread creates a team of threads (fork) to work in parallel. When all the threads end their task, they are joined together. owing to an implicit barrier. This type of construct is also known as non-iterative worksharing construct. Other types of worksharing constructs are available. Although there are implicit communications between threads through access to shared variables or the implicit synchronization at the end of parallel regions and worksharing constructs, explicit synchronization mechanisms for mutual exclusion are also available in OpenMP. These are critical or atomic directives, lock routines, and event synchronization directives. OpenMP is not responsible for the management of the memory hierarchy but certain issues regarding cache memory management should be borne in mind. There are two factors that determine whether a loop schedule is efficient: data locality and workload balancing among iterations. The best schedule that we can choose when there are data locality and a good workload balance is static with a chunk size of q=n/p, where nis the number of iterations and pthe number of threads. In other cases, dynamic or guided schedules may be adequate. When a cache line, shared among different processors, is invalidated as a consequence of different processors writing in different locations of the line, false sharing occurs. False sharing must be avoided as it decreases performance due to cache trashing. One way to avoid this is to divide the data to be accessed by different processors into pieces whose size is multiple of the cache line size. A good practice for improving cache performance is to choose a schedule with a chunk size that minimizes the requests of new chunks and that is a multiple of a cache line size. 40 Chapter 2. Fundamentals 2.8.2 CUDA The NVIDIA Compute Unified Device Architecture (CUDA) is a parallel computing platform and a programming model. The GPU provides massively parallel processing capabilities with a high computational throughput due to their large number of cores. In particular, CUDA is organized into a set of Streaming Multiprocessors (SMs), each one containing many Scalar Processor (SP) with many cores inside. The NVIDIA’s G80 series, introduced in November 2006, had a total of 128 cores in 16 SMs, each one with 8 SPs, as summarized in Figure 2.19(a). This architecture uses a Single Instruction, Multiple Thread (SIMT) parallel programming model. SIMT is the terminology used by NVIDIA to define a hybrid model between vector processing (SIMD) and hardware threading [40]. This approach makes it possible to write single instructions which will be simultaneously executed from multiple threads. The general specification, as well as the number of resources available on the GPU depends on its compute capability (CC), which is represented by a version number, 1.x, 2.x, 3.x and 5.x [119]. The first CUDA capable hardware, codenamed tesla, has a CC of 1.x and defines the main features of the architecture, such as the organization into a set of SMs, SPs and the different types of device memory. A full list of the differences among each compute capability can be found in [119, 40]. The basic compute unit in CUDA is the SM which has fixed and limited resources, such as a set of registers and on-chip memory that are shared among the cores within the same SM. The GPU can manage and schedule thousands of threads in hardware simultaneously, avoiding high thread management overheads. A large number of threads are required to make full use of the GPU computing capabilities. The threads execute the same instruction on different data by grouping 32 threads as the minimum size of collaborative unit, called a warp. The warp size is implementation defined and it is related to shared memory organization, data access patterns and data flow control [90]. The threads are arranged in a 1D, 2D or 3D grid of blocks which are scheduled to any of the available SMs. Figure 2.19(a) shows a 2D grid of 8 blocks each one with 4×4 threads. Each thread has a unique ID that identifies which block the thread belongs to, and the thread’s position within the block. The grid of blocks is usually designed to match a thread with a value among the different data that will be processed on the GPU. One of the most important aspects of the architecture is the memory hierarchy as it plays a key role in performance. As shown in Figure 2.19(a), the G80 architecture (CC 1.0) has a global memory, a texture memory and a constant memory which are available for all the 2.8. Parallel programming models 41 (a) (b) Figure 2.19: CUDA. (a) Overview of G80 architecture, (b) grid of blocks and block of threads scheduled to any of the available SMs. threads at any Streaming Multiprocessor. There was no a L1 / L2 cache hierarchy for caching global memory data on the early GPU architectures. Therefore, the device memory was designed for specific purposes. The global memory is the main device memory on the GPU, and has the slowest access time. However, the number of memory accesses to global memory can be reduced for certain memory access patterns by a technique known as coalescing. If threads within the warp7 request consecutive and aligned values in memory, the device coalesces the global memory transactions (load/store) into as few transactions as possible (one in the best case) [119]. Data allocated in the texture memory space will be automatically cached. There are only 8 KB of cache per SM and it is optimized for 2D spatial locality. This can improve performance when threads access values in some regular spatial neighborhood. Texture memory is generally used for visualization or as input buffers owing to data caching. The constant memory is another possibility of caching accesses (64 KB in total) to global memory. A single read from constant memory can be broadcast to a warp, if the same value is accessed by all the threads within the warp. As a result, reading from this memory costs one memory transaction from global memory the first time data are accessed, and one read from the constant cache otherwise. The on-chip memory, known as shared memory, enables extremely rapid load/store accesses to the data but within the lifetime of the block. There are 16 KB of shared memory 7On hardware of compute capability 1.x, memory transactions are coalesced within half warp. 42 Chapter 2. Fundamentals within each SM but it is only “shared” for the threads of the same block. Therefore, sharing data among different blocks is not possible through this memory and becomes a challenge when programming for the GPU. It is up to the programmer loading data from global (texture) memory to shared memory. The main feature of the shared memory is reusing data within a block, sharing data among the threads of the same block. This way, the shared memory can be managed as an explicit cache defined by the programmer. Although the amount of shared memory per SM is only 16 KB, the effective use of this memory can lead to speedups of 7× compared to a naive implementation in global memory [40]. Each thread on the GPU has its own local memory and a set of registers where the computation takes place. The maximum number of registers that can be used by a thread is limited by the CC of the GPU, as well as the number of threads configured per block [119]. Different aspects must be taken into consideration to efficiently exploit all the memory on the GPU [90]. For example: (1) data transfer between the CPU and the GPU should be minimized; (2) the memory access pattern must be coalesced to consecutive memory locations; (3) data loaded in shared memory should be reused to reduce load/store accesses to the global memory; (4) the number of threads per block must be optimized to run the maximum concurrent threads allowed in each SM; and (5) minimizing the global synchronization among blocks leads to a reduction in the execution time of the program. Paying attention to these aspects is vital for GPU programming. In Fermi (CC 2.x), Kepler (CC 3.x) and Maxwell (CC 5.x) architectures, there is also an L1 and L2 cache hierarchy, as shown in Figure 2.20 for Kepler. The accesses to global memory are cached in this memory hierarchy. The L2 cache is fully available for all the threads and the L1 cache only for the threads running in the same Streaming Multiprocessor (in Kepler, the SMs are called SMXs). Note that the L1 and L2 caches are managed by the GPU, unlike the shared memory, but the programmer can take advantage of this cache hierarchy by exploiting the memory access pattern when reading data from the global memory (the same when writing data to the global memory). The L1 cache is placed in the same chip as the shared memory, so the size of the on-chip memory per SMX (see Figure 2.20) is split between these two kind of memories. In the following, we will cover the major characteristics found in the CC 2.x and 3.x corresponding to the GPUs used in this work. The main differences are summarized in Table 2.1. The main changes introduced in CC 2.x are: (i) extension in size of the shared memory from 16 KB up to 48 KB per SM, (ii) configurable 16 KB or 48 KB of L1 cache on each SM, 2.8. Parallel programming models 43 Figure 2.20: Overview of the Kepler architecture incorporating a L1 and L2 cache hierarchy. and (iii) shared L2 cache for all SMs. The CC 3.x mainly increases the number of resources available per thread, SP and SMX compared to CC 2.x. The CC 3.5 has a read-only cache of 48 KB shared by every three SMXs. In all the architectures, threads within the same block can be synchronized, for example to communicate intermediate results in the shared memory as part of a parallel computation. However, it is not possible to synchronize threads among different blocks. Owing to this restriction, the communication among all the threads must be through the global memory and this becomes another challenge when a thread needs data which have been generated outside its block. 44 Chapter 2. Fundamentals G80 Fermi GF110 Kepler GK104 Kepler GK110 Compute capability 1.0 2.0 3.0 3.5 Max SMs 8 16 8 15 CUDA cores per SM 8 32 192 192* Max threads per block 512 1024 1024 1024 Max total threads per SM 768 1536 2048 2048 Max number of blocks per SM 8 8 16 16 Max Shared Mem 16 kB 48 kB 48 KB 48 KB Max L1 cache N/A 48 kB 48 KB 48 KB Max L2 cache N/A 768 KB 512 KB 1536 KB Max Memory 1536 MB 1536 MB 2048 MB 6144 MB *Plus an additional 64 dual-precision units per SM. Table 2.1: Max GPU resources defined by the compute capability. A CUDA program is called a kernel and it is executed on the GPU by thousands of threads in parallel. A kernel is configured by specifying the grid size, the threads within the block, and the amount of shared memory used by a block. When a kernel is called, the computation on the GPU takes place. The blocks, which are arranged into a grid, are scheduled to any of the available cores enabling automatic scalability for future architectures. 2.9 Experimental setup In this thesis we have taken into account the efficient computation on GPU of the techniques under study. In recent years, GPUs (and the CUDA API) have evolved so rapidly that we have evaluated our work on different architectures (and CUDA versions). In this section we present the main specifications of the CPU and GPUs used in our experiments. The experiments are executed under Linux using the gcc compiler version 4.6.3 for the OpenMP implementations, and the nvcc compiler for the case of the CUDA implementations, respectively, with full optimization flags (-O3) in both cases. We describe the performance measures to evaluate the classification results in terms of accuracy, as well as in terms of execution time and speedup. Different datasets have been used in this work and will be also detailed in this section. The names and number of known samples (reference map) of the hyperspectral images are summarized at the end of this section. 2.9. Experimental setup 45 2.9.1 Hardware used in the experiments In this work we have used an Intel quad-core i7-860 microprocessor (8MB Cache, 2.80 GHz) and 8 GB of RAM as the base architecture for comparison. Each core has a separated L1 cache for instructions and data, and a unified L2 cache. The unified L3 cache is common to all the cores. The main characteristics of the memory hierarchy are described in detail in [82]. Table 2.2 summarizes the main specifications of this CPU. We have used the NVIDIA GTX 580, GTX 680 and GTX TITAN devices. Table 2.3 shows the GPU model, the compute capability of the graphic cards, available resources, and CUDA versions used in this thesis. # cores Core clock RAM L1 L2 L3 (GHz) (GB) (KB) (KB) (KB) Intel core i7-860 4 2.80 8 64/64 256 8192 Table 2.2: Main characteristics of the CPU used in this thesis. GTX 580 GTX 680 GTX TITAN Compute capability 2.0 3.0 3.5 Streaming Multiprocessors 16 8 15 CUDA cores per SM 32 192 192* Threads per block 1024 1024 1024 Total threads per SM 1536 2048 2048 Number of blocks per SM 8 16 16 Active warps per SM 48 64 64 Registers per SM 32768 65536 65536 Registers per thread* 63 63 63 Device Memory 1536 MB 2048 MB 6144 MB Shared Memory 48 KB 48 KB 48 KB L1 cache 48 KB 48 KB 48 KB L2 cache 768 KB 512 KB 1536 KB CUDA version 4.0 5.0 5.5 *Plus an additional dedicated zero register. Table 2.3: GPU model, compute capability, resources and CUDA version of the graphic cards used in this thesis. 52 Chapter 2. Fundamentals Indian Pines Salinas Valley Classes Number Classes Number of samples of samples 1-Alfalfa 54 1-Brocoli_green_weeds_1 2009 2-Corn-notill 1434 2-Brocoli_green_weeds_2 3726 3-Corn-mintill 834 3-Fallow 1976 4-Corn 234 4-Fallow_rough_plow 1394 5-Grass/pasture 497 5-Fallow_smooth 2678 6-Grass-trees 747 6-Stubble 3959 7-Grass/mowed 26 7-Celery 3579 8-Hay-windrowed 489 8-Grapes_untrained 11271 9-Oats 20 9-Soil_vinyard_develop 6203 10-Soybean-notill 968 10-Corn_senesced_green_weeds 3278 11-Soybean-mintill 2468 11-Lettuce_romaine_4wk 1068 12-Soybean-clean 614 12-Lettuce_romaine_5wk 1927 13-Wheat 212 13-Lettuce_romaine_6wk 916 14-Woods 1294 14-Lettuce_romaine_7wk 1070 15-Bld-Grs-Trs-Drs 380 15-Vinyard_untrained 7268 16-Stone-Steel 95 16-Vinyard_vertical_trellis 1807 Total 10366 Total 54129 Table 2.7: Total number of samples for AVIRIS datasets: Indian Pines and Salinas Valley. Hekla Volcano Classes Number Classes Number of samples of samples 1-Andesite lava 1970 342 7-Hyaloclastite formation 684 2-Andesite lava 1980 I 708 8-Lava covered 700 3-Andesite lava 1980 II 1496 9-Rhyolite 404 4-Andesite lava 1991 I 2739 10-Scoria 550 5-Andesite lava 1991 II 410 11-Firn and glacier ice 458 6-Andesite lava with moss 1023 12-Snow 713 Total =10227 Table 2.8: Total number of samples for AVIRIS dataset: Hekla Volcano. The Indian Pines scene was acquired over a mixed agricultural/forested region in northwestern Indiana with a moderate spatial resolution of 20 m. This image represents a very challenging land-cover classification scenario. It consists of 145×145 pixels and 220 spectral 2.9. Experimental setup 53 bands. The four bands covering the region of water absorption were removed. This dataset is available through Purdue’s University MultiSpec site. Figure 2.25(a) shows the true color representation of the scene. The sixteen classes of interest available in the reference map are shown in Figure 2.25(b). The total number of samples are presented in Table 2.7. The number of training samples for this scene is usually randomly taken from the reference map as 5% or 10% of the available data. The Salinas dataset has 512×217 pixels and it was captured over the Salinas Valley in California. It is characterized by a high spatial resolution (3.7 m) owing to a low-altitude flight during the acquisition [72]. None of the hyperspectral bands was removed in this dataset. The total number of samples are presented in Table 2.7, and like the Indian Pines, the training samples are randomly taken from the reference map as 5% or 10% of the available data. Figure 2.26(a) shows the true color representation of the scene, and Figure 2.26(b) shows the reference map. The Hekla volcano scene has spatial dimensions of 560×600 pixels with a spatial resolution of 20 m. During data collection of this dataset, spectrometer four was not properly working. This particular spectrometer operates in the near-infrared wavelength range, from 1.84 µm to 2.4µm (64 data channels). These 64 data channels were deleted from the data set along with the first channels for all the other spectrometers, but those channels were blank. Once the noisy and blank data channels had been removed, 157 data channels were left [16]. Figure 2.27(a) shows a false color representation of the scene. The twelve classes of interest available in the reference map are shown in Figure 2.27(b). The names and the number of samples per class are available in Table 2.7. For this dataset, a fixed number of fifty samples per class is used for training, and the rest of the samples are used for testing. Table 2.9 summarizes the dimensions and the size in MB of the five hyperspectral images used in the experiments carried out in this thesis. We have considered different spatial resolutions (1.3µm, 3.7µm and 20 µm), different dimensions in the spatial domain (from 145×145 to 1096×715 pixels), and different numbers of spectral bands (from 102 to 224 bands). The Pavia City scene is the largest in the spatial domain with 1096×715 pixels, while the scenes of Indian Pines and Salinas Valley are the scenes with more spectral bands, 220 and 224, respectively. The Pavia University, Pavia City, Indian Pines and Salinas Valley datasets (hyperspectral image and reference map) are publicly available online9at the research webpage of the 9http://www.ehu.eus/ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes 54 Chapter 2. Fundamentals Dataset name Dimensions Size (MB) Sensor Pavia University 610×340×103 162.9 ROSIS-03 Pavia City 1096×715×102 609.8 ROSIS-03 Indian Pines 145×145×220 35.3 AVIRIS Salinas Valley 512×217×224 189.9 AVIRIS Hekla Volcano 560×600×157 402.5 AVIRIS Table 2.9: Hyperspectral images used in this work. Computational Intelligence Group from the Basque University (UPV/EHU). The Indian Pines scene is originally available through Purdue’s University MultiSpec site10. We would like to thank prof. Benediktsson from the University of Iceland, for providing the hyperspectral dataset of Hekla. 10http://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral.html 2.9. Experimental setup 55 (a) (b) (c) Figure 2.23: Pavia University dataset. True color representation (a), reference map with nine classes of interest (b), and name of the classes (c). 56 Chapter 2. Fundamentals (a) (b) (c) Figure 2.24: Pavia City dataset. True color representation (a), reference map with nine classes of interest (b), and name of the classes (c). 2.9. Experimental setup 57 (a) (b) (c) Figure 2.25: Indian Pines dataset. True color representation (a), reference map with sixteen classes of interest (b), and name of the classes (c). (a) (b) (c) Figure 2.26: Salinas Valley dataset. True color representation (a), reference map with sixteen classes of interest (b), and name of the classes (c). 58 Chapter 2. Fundamentals (a) (b) (c) Figure 2.27: Hekla Volcano dataset. False color representation (a), reference map with twelve classes of interest (b), and name of the classes (c). CHAPTER 3 SPECTRAL-SPATIAL CLASSIFICATION SCHEMES BASED ON SEGMENTATION AND MATHEMATICAL MORPHOLOGY 3.1 Introduction In this chapter we present the two schemes we proposed for efficient spectral-spatial nD image classification, based on segmentation, Mathematical Morphology (MM), and SVM classifiers, called CA–WSHED–MV and WT–EMP. First, we describe in Section 3.1.1 the general framework for spectral (pixel-wise) classification schemes as a basis for designing new schemes. This framework is a starting point for hyperspectral image classification. Second, a general scheme for spectral-spatial classification, as well as the common data fusion strategies for joining the spectral and the spatial information of such scheme are described in this chapter in Section 3.1.2. This framework incorporates spatial information to improve the results of the final classification. Based on the approach used for extracting the spatial information, different data fusion techniques can be employed. The spectral-spatial CA–WSHED–MV scheme, originally proposed in [155], is presented in Section 3.2. In this scheme extracts the spatial information by a watershed transform based on cellular automata (CA–Watershed), resulting in a more efficient step for GPU processing, combining the spectral results of the classifier by a Majority Vote (MV). 60 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM The spectral-spatial WT–EMP scheme is a new proposal for land-cover classification in remote sensing imaging for real-time applications. The second scheme proposed in this thesis extracts spatial information using morphological profiles and creates a new vector of features prior to the classification. The scheme is described in Section 3.3. The classification accuracy of both schemes is evaluated on real hyperspectral images, and the results thereof are compared to classification schemes found in the literature based on segmentation and MM. 3.1.1 Framework for spectral classification schemes The classification schemes for hyperspectral images process each pixel independently (pixelwise processing) and do not take into account the spatial information of the neighborhood. Success in classification depends solely on the ability of the classifier to discriminate among pixel vectors. Thus, the classifier discriminates the pixels of the image based only on the spectral features captured by the hyperspectral sensor. One possible framework for spectral classification that we will use in this thesis is given in Figure 3.1. This framework is based on the spectral classification scheme proposed by Landgrebe et al. [95], which is widely used nowadays as the basis for hyperspectral image classification. Considering that the spectral features are often redundant, Feature Extraction (FE) and Feature Selection (FS) are usually performed as a preprocessing step. FE is designed to remove the redundancy introduced by the spectral correlation between bands [146], while the main characteristics in the spectral domain are retained. One consequence is that the spectral dimensionality is reduced. Different techniques have been investigated for extracting the principal features for urban land cover classification of hyperspectral data using SVM-based classifiers [48, 31]. Section 2.2.1 describes the techniques used in the spectral-spatial classification schemes investigated in this work. Principal Component Analysis (PCA) and Independent Component Analysis (ICA) are two of the most widely used unsupervised techniques in remote sensing for reducing the redundancy of the spectral bands. However, the Minimum Noise Fraction (MNF), which extracts the data sorting the components by signal-to-noise ratio, performed better than PCA in the analysis carried out in [48]. This finding emphasizes the idea that removing noise is a good preprocessing option. Supervised feature extraction techniques such as Discriminant Analysis Feature Extraction (DAFE), Decision Boundary Feature Extraction (DBFE) and Non-parametric Weighted Feature Extraction (NWFE) have shown their effectiveness in classification of hyperspectral 3.1. Introduction 61 Figure 3.1: Simplified framework for spectral classification schemes incorporating feature reduction and denoising. data [108, 31, 93, 96]. Kuo and Landgrebe proposed NWFE in [93], as a method for improving DAFE by exploiting different weights on samples close to the decision boundary. The NWFE has been shown to perform better than DAFE and DBFE in classification of hyperspectral data [96]. An automatic extraction of features using wavelet was presented in [86], where the dimensionality of the image is reduced several times, and then reconstructed by wavelets. The best level of decomposition is then determined through the similarity between the original pixel vector and the vector reconstructed using wavelets. In [162] several methods based on wavelet transforms are developed to extract useful features for classification. The Daubechies 3 wavelet is applied to hyperspectral images, and a small subset of wavelet coefficients is then used to extract the effective features for classification. In the work published in [68] wavelet transforms were used to reduce the dimension of the hyperspectral images, owing to its ability to detect the local energy variations better than other transforms. The acquisition of the image usually introduces undesirable artifacts that may affect the quality of the spectral signature. Therefore, another preprocessing commonly applied are denoising [164, 166] and scatter correction [140]. As illustrated in Figure 3.1, feature extraction and denoising are not mutually exclusive. Finally, once the preprocessing has been performed, the classification of the hyperspectral image takes place. The result is a thematic map in which each pixel is labelled with a class that identifies the corresponding spectral signature. In Figure 3.1, the result is a classification map where the labels are identified by three different colors. There are a variety of classifiers, such as maximum likelihood, neural networks, decision trees, and kernel-based methods. Melgani and Bruzzone [110] have shown that the SVM classifier is more effective for remote sensing classification than other conventional non-parametric classifiers, in terms of accuracy, computational time, and stability to parameter setting. In addition, the SVM classifier has been 68 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM 8-connectivity with s=1. The result is a one-band gradient image which will be segmented by the CA–Watershed algorithm. 3.2.2 Watershed transform based on cellular automata (CA-Watershed) Regarding segmentation, Galilée et al. [64] introduced a parallel algorithm-architecture based on asynchronous CA to compute the watershed transform. The asynchronous behavior of this automaton is subject to its implementation; i.e., the CA–Watershed can be synchronously or asynchronously implemented. The 3–state cellular automaton proposed in [64] follows the Hill-Climbing algorithm described in Section 2.3.2. The main advantage of this algorithm is that minima detection, labelling, and climbing the steepest paths are performed simultaneously and locally. We now go on to describe the CA–Watershed algorithm, where each cell of the automaton computes a pixel of the image. The fundamentals of the watershed transform and cellular automata are described in Section 2.3.2 and Section 2.4, respectively. Figure 3.6 shows the 3–state cellular automaton that computes the watershed transform which is described hereafter. First, the pixels are sequentially labelled. In the initial state, all pixels compute the sets NF(u)and N=(u)corresponding to an arbitrary lower slope and the neighbors with the same grey value, respectively. If a pixel has no lower slope, NF(u) = /0, it switches to the minimum or plateau (MP) state and its label is initialized as follows: l(u) = min(l(v)),(3.1) Compute and non-minimum MP NM INIT look at look at Extend plateaus Hill-Climbing Figure 3.6: 3-state cellular automaton implementing Hill-Climbing algorithm [64]. 3.2. Scheme based on segmentation: CA–WSHED–MV 69 Algorithm 1 CA–Watershed algorithm [64] 1: procedure CA–WATERSHED(state,f) 2: switch state do 3: case INIT : .Initialize automaton 4: compute N=(u),NF(u) 5: if (NF(u) = /0) then 6: l(u)←minv∈N=(u)(l(v)) .Eq. (3.1) 7: state(u)←“MP” 8: else 9: f(u)←f(NF(u)) .Eq. (3.2) 10: state(u)←“NM” 11: end if 12: 13: case MP : .Minimum or Plateau 14: for each pixel v∈N=(u)do 15: if h(v)<h(u)then 16: NF(u) = {v}.Eq. (3.4) 17: f(u)←f(v) 18: state(u)←“NM” 19: else 20: l(u)←min(l(u),l(v)) .Eq. (3.3) 21: end if 22: end for 23: 24: case NM : .Non Minimum 25: v=NF(u) 26: f(u)←f(v).Eq. (3.5) 27: 28: end switch 29: end procedure where l(v)are the labels of the pixels v∈N=(u). Otherwise, the state of the pixel switches to non-minimum (NM) and NF(u)points to the climbing direction. The grey value and the label of the pixel are modified as follows: f(u) = fNF(u),(3.2) where f(u)is the pair h(u),l(u),h(u)is the grey value of the pixel uand l(u)its label. Once the pixel has been initialized, the update stage begins. This is an iterative task that processes the MP and NM states. A pixel in MP state waits for data from any neighbor v∈N=(u), and, depending on the grey value of the neighbor, two cases are considered. If the value of the neighbor is equal to or greater than its current value, the label of the pixel is 70 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM updated as: l(u) = min(l(u),l(v)) ifh(v)≥h(u).(3.3) Otherwise, the pixel belongs to the lower border of the plateau and its state switches to NM as follows: NF(u) = v, f(u) = f(v), state =NM,        ifh(v)<h(u).(3.4) On the other hand, a pixel in NM state remains in that state and waits for data from the neighbor NF(u), updating its data as follows: f(u) = f(NF(u)).(3.5) This algorithm is non-deterministic and may lead to different segmentation results. A formal proof of correctness and convergence towards a watershed segmentation using a mathematical model of data propagation in a graph is presented in [64]. The pseudocode of this algorithm is listed in Algorithm 1 and Figure 3.7 shows a guided example of the evolution (a) (b) (c) (d) (e) (f) Figure 3.7: CA–Watershed based on Hill-Climbing algorithm. Example of 1D image represented as a terrain with dashed lines and as grey values with the squares: a) pixels are sequentially labelled and states are initialized, b) MP updating stage where two pixels change their state to NM c) to f) the automaton evolves only with the NM rules. 3.2. Scheme based on segmentation: CA–WSHED–MV 71 of the CA–Watershed on a 1D image, representing the grey values with the squares and the terrain with dashed lines. The automaton initialization is shown in Figure 3.7(a). The first updating step, Figure 3.7(b) and Figure 3.7(c), illustrates the MP and NM states, respectively. It can be observed in Figure 3.7(b) that two pixels change their state from MP to NM as they belong to the lower border of the inner plateau. In the successive updating steps, Figure 3.7(d) to Figure 3.7(f) the automaton evolves only with the NM rules. The labels are propagated from the minimum climbing up the hills, flooding the terrain. 3.2.3 Results In this section we conduct an analysis1. of the proposed CA–WSHED–MV scheme in terms of Overall Accuracy (OA), Average Accuracy (AA) and the kappa coefficient of agreement (k). These criteria are computed from the confusion matrix. The RCMG is computed using 8-connectivity and the pairs of pixel vectors removed is s=1, (see Section 2.2.2 and Section 3.2.1 for details on the robustness of this morphological gradient). The CA–Watershed is configured to use 8-connectivity as well. Given that the classification results depend on the set of training samples and their distribution within the hyperspectral image [39], the results are obtained executing the classification 10 times with different sets of training samples each time. Since the training samples were randomly selected, the results were calculated as the average of the 10 different values obtained for each execution. The hyperspectral remote sensing scenes used in the experiments are two urban areas (Pavia University and Pavia City datasets) taken by the ROSIS-03 hyperspectral sensor and one hyperspectral image of a crop area (Indian Pines dataset) taken by the AVIRIS sensor. These datasets are described in Section 2.9.3, including the number of available samples. The number of training samples for each scene is chosen as in [155, 121, 54] to compare our scheme on equal terms. In order to obtain a reliable evaluation of the results, the accuracies are calculated excluding the samples used for training. 1Part of these results have been published in P. Quesada-Barriuso, F. Argüello, and D. B. Heras, “Efficient segmentation of hyperspectral images on commodity GPUs,” in 16th International Conference on Knowledge-Based and Intelligent Information & Engineering System, vol. 243, pp. 2130–2139, 2012. And P. Quesada-Barriuso, F. Argüello, and D. B. Heras, “Computing efficiently spectral-spatial classification of hyperspectral images on commodity gpus,” in Recent Advances in Knowledge-based Paradigms and Applications (J. W. Tweedale and L. C. Jain, eds.), vol. 234 of Advances in Intelligent Systems and Computing, Ch. 2, pp. 19–42, Springer International Publishing, 2014. 72 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM Pavia University Pavia City Indian Pines C128 128 1024 γ0.125 0.125 0.5 Table 3.1: Best parameters (C,γ)for the SVM determined for each hyperspectral image in the range C=[1, 4, 16, 64, 128, 512, 1024], γ=[0.5, 0.25, 0.125, 0.0625] by 5–fold cross-validation. The whole process is computed in CPU and the classification is carried out using the LIBSVM library [35] with the RBF kernel. The parameters (C,γ)for training the SVM are determined for each hyperspectral image in the range C=[1, 4, 16, 64, 128, 512, 1024], γ= [0.5, 0.25, 0.125, 0.0625] by 5–fold cross-validation, as explained in Section 2.7. For each one of the hyperspectral images under study, a random training set of the available samples is used to find the best parameters for the SVM, which are summarized in Table 3.1. The classification results are compared to those obtained by the pixel-wise SVM classifier and other spectral-spatial classification schemes based on segmentation and SVM, the results of which are available in the literature: SVM+EM (PR) [154], SVM+ISODATA [154], SVM+ ISODATA (PR) [154], and HSeg+MV [158]. These schemes are based on clustering and region growing segmentation techniques. The HSeg+MV scheme is based on the HSEG algorithm proposed in [161]. A description of these segmentation techniques was presented in Section 2.3. The schemes marked with PR also include an additional post regularization as described in the previous Section. The segmentation map obtained by these schemes is combined with the thematic map produced by the SVM using the majority voting data fusion. These segmentation techniques are also described in Section 2.3. For the Pavia City dataset there are no published results of spectral-spatial classification schemes based on segmentation2, so we have produced additional results for this dataset by using two segmentation algorithms based on clustering. The first classification scheme is based on k-means (k-means–MV) and the second one is based on Quick-Shift [163] (QS– MV). Both schemes combine the spectral and spatial information by majority vote. University of Pavia dataset The University of Pavia dataset is a moderately dense urban area, with some buildings and large meadows and high spatial resolution (1.3 m). Nine classes of interest are available for 2Partial result for this dataset were found in [21] but using a subset of 900×300 pixels. 3.2. Scheme based on segmentation: CA–WSHED–MV 73 SVM SVM+EM SVM+ISODATA Hseg+MV CA–WSHED–MV [158] (PR) [154] (PR) [154] [158] 1-Asphalt 84.93 93.45 94.40 94.77 95.42 2-Meadows 79.79 93.78 87.45 89.32 95.51 3-Gravel 67.16 82.53 61.32 96.14 87.01 4-Trees 97.77 99.38 98.63 98.08 95.93 5-Metal sheets 99.46 100 99.91 99.82 99.72 6-Bare Soil 92.83 97.38 97.88 99.76 97.26 7-Bitumen 90.42 94.19 100 100 98.97 8-Bricks 92.78 98.31 99.02 99.29 98.51 9-Shadows 98.11 97.86 97.86 96.48 100 OA 81.01 94.64 91.20 93.85 95.88 AA 88.25 95.21 92.94 97.07 96.48 k0.758 0.929 0.884 0.919 0.944 Table 3.2: Classification results for the CA–WSHED–MV scheme on the University of Pavia dataset and compared to SVM, SVM+EM (PR), SVM+ISODATA (PR), and Hseg+MV. Best results are indicated in bold. this scene. The number of training samples used in the experiments was taken from [54], as described in Section 2.9.3, Table 2.6. We recall that the number of training samples is not included in the test set. Table 3.2 shows the name of the classes and gives the classification accuracies obtained by the pixel-wise SVM classifier, and different segmentation schemes based on clustering and post-regularization SVM+EM (PR) and SVM+ISODATA (PR), and based on Hierarchical Image Segmentation Hseg+MV. Best results are indicated in bold. Figure 3.8 shows from left to right, the thematic map produced by the pixel-wise SVM classifier, the RCMG, the CA– Watershed segmentation map with watershed lines imposed for visualization, and the final classification map produced for the CA–WSHED–MV scheme. As can be seen from Table 3.2, the spectral-spatial classification schemes have higher accuracies as compared to the results obtained by the pixel-wise SVM. The best OA is achieved by the proposed CA–WSHED–MV scheme. This scheme improves the OA in 14.87% and the AA in 8.23% as compared to the classification by the pixel-wise SVM. The Hseg+MV scheme has the best AA, showing a good consistency in the class-specific results, which were improved for almost all the classes, except for the shadows class. However, the proposed scheme improved all the class-specific accuracies, including the shadows class up to 100%. This is an outstanding result because the accuracies in the sixth column in Table 3.2, corresponding to the proposed scheme, are calculated as the average over 10 different classification with different sets of training and test samples. 74 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM Pavia City dataset The Pavia City dataset is a dense urban area, with many buildings and a large river in the top-right corner of the scene. The nine classes of interest, which are available for this dataset are shown in Table 3.3. The number of samples used in the experiments is detailed in Section 2.9.3, Table 2.6. Table 3.3 shows the classification accuracies obtained by the pixel-wise SVM classifier. Other schemes based in segmentation applied to this dataset in the same experimental conditions were not found in the literature. K. Bernard et al. [21] published the results of their spectral-spatial classification scheme based on a stochastic minimum spanning forest approach. Their scheme produces a set of segmentation maps generating random markers from the SVM thematic map, which are aggregated by a MV. However, the Pavia City dataset used in [21] is a subset of 900×300 pixels, so it is not included in the analysis because the number of test samples would not be the same. In order to obtain a reliable evaluation of the proposed scheme, we have added the results of a partitional clustering approach by k-means (k-means–MV) generated by us. We have also used a Quick Shift (QS) segmentation technique [163] to create a second clustering-based scheme (QS–MV) for comparison. These schemes are similar to the one proposed by Tarabalka et al. [154] but without post-regularization. Figure 3.9 shows from left to right the thematic map produced by the pixel-wise SVM classifier, the RCMG, the CA–Watershed segmentation map with watershed lines imposed for visualization, and the final classification map produced for the proposed CA–WSHED–MV scheme. SVM k-means–MV QS–MV CA–WSHED–MV [58] 1-Water 99.1 99.99 99.97 100 2-Trees 90.8 92.84 96.29 95.84 3-Meadow 97.4 63.25 96.89 98.33 4-Brick 87.5 97.49 98.65 98.47 5-Bare Soil 94.6 94.49 94.79 94.76 6-Asphalt 96.4 97.85 94.59 99.09 7-Bitumen 96.5 92.65 96.69 95.22 8-Tile 99.5 98.88 99.36 99.80 9-Shadows 100 95.53 94.58 100 OA 98.1 97.93 98.76 99.20 AA 95.1 92.56 96.84 97.94 k0.97 0.970 0.982 0.988 Table 3.3: Classification results for the CA–WSHED–MV scheme on the Pavia City dataset and compared to SVM, k-means–MV and QS–MV. Best results are indicated in bold. 3.2. Scheme based on segmentation: CA–WSHED–MV 75 Table 3.3 gives the classification accuracies obtained by the pixel-wise SVM classifier, and the segmentation schemes based on k-means and QS, as well as the proposed scheme. Results for the pixel-wise SVM classifier are taken from [58]. As can be seen from Table 3.3, all the spectral-spatial classification accuracies are higher as compared to the SVM pixel-wise classifier, despite the OA obtained as basis for comparison is fairly high (98.1%). The best OA is achieved when using the spectral-spatial classification scheme based on segmentation by the CA–Watershed algorithm. With the proposed scheme, the OA is improved by 1.1% and the AA is improved by 2.87% compared to the pixel-wise SVM classification. The kmeans–MV scheme obtains worse results mainly because the class-specific accuracy for the class meadows is only of 63.25, which is 35.15% less than the base for comparison. The QS–MV shows a good regularity in the class specific accuracies with an AA of 96.84% but the CA–WSHED–MV scheme is more regular among the different classes (AA =97.94%). Indian Pines dataset The Indian Pines dataset has a low spatial resolution (20 m/pixel) and similar agricultural and forested regions. Sixteen classes of interest are available in the reference map. The total number of samples is presented in Section 2.9.3, Table 2.7. The number of samples used for training the SVM is randomly taken from the reference map as 10% of the available data as in [154, 58]. Table 3.4 gives the classification accuracies obtained by the pixel-wise SVM classifier, the segmentation scheme based on ISODATA [154] and the segmentation scheme based on hierarchical image segmentation (HSeg+MV) [58]. Best results are indicated in bold. Figure 3.8 shows from left to right the thematic map produced by the pixel-wise SVM classifier, the RCMG, the CA–Watershed segmentation map with watershed lines imposed for visualization, and the final classification map produced for the CA–WSHED–MV scheme. The Hseg+MV scheme obtained the best OA (90.8%), AA (94.0%) and kcoefficient (0.90%) for this scene, as shown in Table 3.4. The proposed scheme has the next highest AA (91.20%) showing also good regularity in the class specific results. The SVM+ISODATA scheme (with and without PR) shows a less uniform result in the CS accuracies, leading to a lower AA as compared to the Hseg+MV and CA–WSHED–MV schemes. This is because the alfalfa class and the oats class are poorly classified by the SVM+ISODATA scheme. Moreover, the post-regularization (fourth column in Table 3.4) reduces the accuracy of these two classes even further. 76 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM SVM SVM+ISODATA SVM+ISODATA Hseg+MV CA–WSHED–MV [154] [154] (PR) [154] [58] 1-Alfalfa 32.65 12.24 6.12 92.3 95.83 2-Corn-notill 74.59 79.32 80.48 90.5 85.22 3-Corn-mintill 64.58 84.95 88.02 83.0 86.14 4-Corn 58.77 75.83 85.31 95.7 93.90 5-Grass/pasture 87.05 93.08 93.75 94.4 96.26 6-Grass-trees 92.72 94.80 96.88 97.6 97.27 7-Grass/mowed 29.17 91.67 91.67 100 100 8-Hay-windrowed 96.37 97.51 97.51 99.5 99.38 9-Oats 22.22 16.67 11.11 100 65.55 10-Soybean-notill 69.76 83.85 84.19 92.1 90.43 11-Soybean-mintill 79.21 93.16 95.77 84.1 75.28 12-Soybean-clean 75.41 85.17 89.87 95.4 95.19 13-Wheat 90.58 93.19 98.43 98.2 99.57 14-Woods 91.07 97.17 97.85 98.6 93.23 15-Bld-Grs-Trs-Drs 65.50 79.53 85.38 82.1 92.28 16-Stone-Steel 84.88 86.05 87.21 100 93.64 OA 78.76 88.53 90.64 90.8 87.95 AA 69.66 79.01 80.60 94.0 91.20 k0.757 0.869 0.893 0.90 0.864 Table 3.4: Classification results for the CA–WSHED–MV scheme on the Indian Pines dataset and compared to SVM, SVM+ISODATA, SVM+ISODATA (PR), and Hseg+MV. Best results are indicated in bold. 3.2.4 Final discussion In this section we have described and analyzed the proposed spectral-spatial n-dimensional image classification scheme based on segmentation and Majority Vote (MV): CA–WSHED– MV. The scheme reduces the dimensionality of the hyperspectral image computing a Robust Color Morphological Gradient (RCMG), and then incorporates the spatial information by a watershed transform based on cellular automata (CA-Watershed). The classification was carried out by SVM combining the spectral and spatial results with a MV between the thematic map produced by the SVM and the segmentation map created by the CA–Watershed. This scheme was originally proposed in [155] but the watershed algorithm used in our scheme has the advantage that it does not create the watershed lines, eliminating the need of a postprocessing of the segmentation map to assign a region to the pixel in the watershed lines. The classification accuracy of the CA–WSHED–MV scheme was proven, producing good results on two urban areas acquired by the ROSIS-03 sensor (University of Pavia and Pavia City datasets), and one mixed agricultural and forested scene (Indian Pines) from the AVIRIS sensor. The classification accuracies were compared to spectral-spatial classification schemes 3.2. Scheme based on segmentation: CA–WSHED–MV 77 found in the literature based on the same framework (segmentation + majority vote): SVM+EM, SVM+EM (PR), SVM+ISODATA, SVM+ISODATA (PR) and Hseg+MV. In addition, other segmentation techniques such as k-means and quick-shift (QS) were used to obtain a reliable evaluation of the proposed scheme in the cases where the same spectral-spatial framework was not found in the literature for comparison. Experimental results have shown that the proposed scheme based on CA–Watershed improves the classification accuracies and performs better than other segmentation-based schemes in urban areas. By including spatial information from the closest neighbors, through a fixed or an adaptive post-regularization, small spatial structures may disappear being assimilated by larger structures [155]. However, the CA–WSHED–MV scheme has shown to be robust to this drawback, showing a good consistency in the class-specific accuracies. In the thematic maps produced by the proposed scheme, see Figure 3.8(d), Figure 3.9(d), and Figure 3.10(d), the spatial structures are more homogeneous as compared to the SVM pixel-wise thematic map. 84 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM h g1g2 0.14301535070442 -0.01850334430500 -0.04603639605741 0.51743439976158 -0.06694572860103 -0.16656124565526 0.63958409200212 -0.07389654873135 0.00312998080994 0.24429938448107 0.00042268944277 0.67756935957555 -0.07549266151999 0.58114390323763 -0.46810169867282 -0.05462700305610 -0.42222097104302 0 Table 3.6: Set of filters used by the 2D-DWT for denoising [148]. 3.3.4 Results In this section the WT–EMP scheme proposed for efficient spectral-spatial classification of hyperspectral images is evaluated in terms of Overall Accuracy (OA), Average Accuracy (AA) and kcoefficient, which are computed from the confusion matrix3. The scheme is created by applying the Cohen-Daubechies-Feauveau 9/7 wavelet (CDF97) wavelet decomposition several times in the spectral domain. The number of levels of decomposition is performed to a fixed level to produce always 4-band coefficients. From the remaining 4-band coefficients, the EMP is created applying 4 openings and 4 closings by reconstruction, using a disk of radius of increasing size r∈[1,3,5,7]. Thus, the size of the EMP is 36 (2n+1×4 band-coefficients), independently of the size of the hyperspectral data. See Section 3.3.2 for more details. For image denoising, four levels of the Double-Density DWT are applied to the hyperspectral data as described in Section 3.3.3. The best threshold is found experimentally for each dataset. The denoised data and the EMP are combined via stack vectors (see Section 3.1.3). The whole process is computed in MATLAB4. The hyperspectral remote sensing scenes used in the experiments are two urban areas (Pavia University and Pavia City datasets) taken by the ROSIS-03 hyperspectral sensor and one hyperspectral image of a crop area (Indian Pines dataset) taken by the AVIRIS sensor. The number of training samples for each scene are selected as in [54, 41] to compare our scheme on equal terms (see Section 2.9.3, Table 2.6 and 2.7). In order to obtain a reliable evaluation of the results, the accuracies are calculated excluding the samples used for training. 3Part of these results have been published in P. Quesada-Barriuso, F. Argüello, and D. B. Heras, “Spectralspatial classification of hyperspectral images using wavelets and extended morphological profiles,” Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, vol. 7, no. 4, pp. 1177–1185, 2014. 4The Double-Density DWT software used in these experiments is available at [147]. 3.3. Scheme based on morphological profiles: WT–EMP 85 Pavia University Pavia City Indian Pines C128 128 1024 γ0.125 0.125 0.5 Table 3.7: Best parameters (C,γ)for the SVM determined for each hyperspectral image in the range C=[1, 4, 16, 64, 128, 512, 1024], γ=[0.5, 0.25, 0.125, 0.0625] by 5–fold cross-validation. In this analysis the results are obtained executing the classification 100 times with different sets of training samples each time. Since the training samples were randomly selected, the results were calculated as the average of the 100 different values obtained for each execution. The classification is carried out using the LIBSVM library [35] with the RBF kernel. The parameters (C,γ)for training the SVM are determined by 5–fold cross-validation and are shown in Table 3.7. We also conduct a standalone analysis of the different stages of the proposed scheme to ascertain whether the EMP created using wavelets and the contribution of the denoising as a preprocessing stage are efficient. The stages are named: 1. WT–EMP†is the stage considering only the EMP created from wavelet coefficients as input to the SVM classifier, as indicated in the upper branch in Figure 3.11. 2. WT–EMP‡is the stage considering only the Double-Density DWT applied to the hyperspectral data, see bottom branch in Figure 3.11. In addition, we have evaluated a variant of the proposed scheme, hereinafter WT–EMP?, by creating the stack vector from the original data without denoising and the EMP created from wavelet coefficients. It is similar to the scheme proposed by Fauvel et al. [54]. The different stages of the proposed scheme are also evaluated separately in presence of noise. The classification results are compared to those obtained by the pixel-wise SVM classifier and other spectral-spatial classification schemes, which results are available in the literature for the datasets used in these experiments: EMP-PCA [58], EMP-ICA [121], Spec-EMP [54], WSHED+MV [58], HSeg+MV [58], DB-DB [108], EMD-DWT [68], KPCAp[108] and NWNW [108]. A brief description of these schemes is given in the following sections where the mentioned schemes are included for comparison. 86 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM Pavia University dataset The University of Pavia is a moderately dense urban area with nine classes of interest available in the reference map. This dataset has 103 spectral bands. For this scene, the spectral dimensionality is reduced 6 times by the 1D-DWT, producing a EMP of 36 features. It has been experimentally found that the best parameters for denoising this hyperspectral image by the Double-Density DWT [148] was by using a threshold of 4. The new pixel vector is stacked, resulting in 139 features which are used to classify the image by the SVM. Table 3.8 gives the classification accuracies obtained by the pixel-wise SVM classifier, EMP-PCA, Spec-EMP, DB-DB, EMD-DWT, WT–EMP†stage, WT–EMP?, as well as the proposed WT–EMP scheme. In the following, these schemes are briefly described: – EMP-PCA [58] is a spectral-spatial classification scheme based on EMPs where the morphological profiles are created from the first principal components extracted by PCA. – Spec-EMP [54] is created as the EMP-PCA scheme and stacks the original spectral information within the EMP. – DB-DB [108] is a spectral-spatial classification scheme based on EMAPs where each attribute profile is created from the first components extracted by DBFE. The attributes considered in this scheme are the area and the standard deviation. The EAPacreated from the area and the EAPscreated from the standard deviation are concatenated via stack vectors to create the EMAP, which is reduced a second time by DBFE. – EMD-DWT [68] uses spatial and spectral features combining Empirical Mode Decomposition (EMD) with wavelets, in order to generate the smallest set of features that leads to the best classification accuracy. Thus, data fusion is implicitly performed by a sequence of techniques applied in the spatial (EMD) and the spectral (DWT) domain. The best results in Table 3.8 are indicated in bold. All results are obtained with the number of training and test samples described in Section 2.9.3, except for the EMD-DWT scheme that considers 1% more of training samples, as reported in their work. By considering only the EMP created from the wavelet coefficients (fourth column in Table 3.8) the improvement of the OA is 12.89%, reaching 93.9%, as compared to the pixel-wise SVM. The WT–EMP†stage as compared to the EMP-PCA improves the OA in 14.8%. All classes except the bare soil class are improved in this stage when compared to the SVM classifier, unlike the EMP-PCA, where neither the gravel class and shadows class are improved (see Table 3.8). The EMP from 1D-DWT (WT-EMP†) performs better than the EMP from PCA (EMP-PCA) in terms of classification accuracy for this dataset. 3.3. Scheme based on morphological profiles: WT–EMP 87 SVM EMP-PCA WT-EMP†DB-DB spec-EMP WT-EMP?EMD-DWT WT-EMP [54] [58] [114] [54] [68] 1-Asphalt 84.5 94.5 97.6 95.43 95.3 98.0 — 98.8 2-Meadows 66.2 72.8 91.6 95.88 73.5 97.2 — 98.8 3-Gravel 72.0 53.2 95.7 100 65.9 96.2 — 98.8 4-Trees 98.0 98.9 98.9 90.94 99.2 99.0 — 99.2 5-Metal sheets 99.5 99.5 99.7 100 99.5 99.7 — 99.9 6-Bare Soil 93.1 58.1 89.9 98.41 84.1 96.3 — 98.5 7-Bitumen 91.2 96.1 96.9 99.18 97.2 97.5 — 99.0 8-Bricks 92.3 95.3 96.1 98.65 96.1 96.7 — 98.0 9-Shadows 96.6 91.2 99.9 99.96 93.5 99.9 —99.9 OA 79.5 79.1 93.9 97.89 83.5 97.4 99.0 98.8 AA 88.1 84.3 96.3 97.60 89.4 97.8 — 99.0 κ0.74 0.73 0.91 0.972 0.79 0.96 0.97 0.98 Table 3.8: Classification results for the WT–EMP scheme on the University of Pavia scene compared to SVM, EMP-PCA, EMAP, spec-EMP, EMD-DWT. Note that in EMD-DWT 1% more training samples are used in the classification. The best accuracies are indicated in bold. The †indicates that the EMP is used as the only input to the SVM and ?that the stacked vector is built from the original data and the EMP. 88 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM (a) (b) (c) Figure 3.12: WT–EMP on the University of Pavia dataset: (a) true color representation, (b) SVM classification map, and (c) WT–EMP classification map. The results of the DB-DB, which uses different attributes to create an extended multiattribute profile, are higher, outperforming the aforementioned schemes. Two classes, gravel and metal sheets, have a class-specific accuracy of 100%. However, the accuracy of the trees class falls 7% as compared to the pixel-wise SVM classifier. By incorporating the original spectral data into the EMP via stack vectors (seventh column in the table), all the class-specific accuracies are improved. This is the same approach as the one proposed in [54] (Spec-EMP) but using 1D-DWT instead of PCA to compute the EMP. The WT-EMP?improves the OA 3.5% and the AA 1.5% as compared to the WTEMP†stage. These results are better not only over the image as a whole, but also over each class-specific accuracy. This indicates that the classification using the EMP created from the features extracted by wavelets performs better when the spectral information is incorporated to the scheme. The OA of 97.4% obtained by the WT-EMP?approach is high, but there is room for improvement by applying the 2D-DWT denoising to each hyperspectral band. If the original spectral data if denoised before combining the results with the EMP (ninth column in Ta- 3.3. Scheme based on morphological profiles: WT–EMP 89 ble 3.8), the best result in terms of class-specific accuracy, AA (99.0%) and kcoefficient (0.98) are obtained. These result correspond to the proposed scheme. In the WT–EMP scheme, all the class-specific accuracies are among the highest, particularly for the asphalt, meadows,trees and bare soil classes, where some of the other schemes presented a number of difficulties in improving those classes. Compared to the EMD-DWT scheme, the OA values are close in both schemes (0.2% higher for the EMD-DWT scheme) with the advantage that the WT-EMP is computationally less costly. Figure 3.12 shows from left to right, the true color representation, the SVM classification map, and the WT–EMP classification map. Pavia City dataset The Pavia City dataset is a dense urban area with 102 spectral bands and nine classes of interest. For this image, a threshold of 4 for denoising the hyperspectral is used, and the spectral dimensionality is reduced 5 times by the 1D-DWT. Thus, the size of the pixel vectors is 138 (123 spectral bands + 36 features from the EMP). The number of training samples used in the experiments is described in Section 2.9.3, Table 2.6. The classification accuracies obtained by the propose scheme, the EMP-ICA [121], the Spec-EMP [54], as well as WT–EMP†stage and the WT–EMP?approach are given in Table 3.9. The accuracy of the pixel-wise SVM is included as a basis for comparison. The EMP-ICA is a spectral-spatial classification scheme based on EMPs where the morphological profiles are created from the first principal components extracted by ICA. In the Spec-EMP scheme, the morphological profiles are created from PCA and then they are stacked to the original spectral information. It can be observed in the table that the base classification accuracy is already fairly high (97.7%) using only the spectral information. The most significant improvements are achieved by the WT–EMP scheme. Comparing the third and fourth columns (two similar approaches), 7 from the 9 classes are better classified with the WT–EMP†stage as compared to the EMP-ICA scheme. The same happens when analyzing the results obtained by the WT–EMP?approach, where all the class specific accuracies are improved. This indicates that the EMP created from the wavelet coefficients is working properly and it is the factor that contributes most to the scheme. The best OA (99.7%), AA (99.3%), and kcoefficient (0.99) are obtained for the proposed scheme, when the denoised data are stacked with the EMP created from wavelets. Figure 90 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM SVM EMP-ICA WT-EMP†spec-EMP WT-EMP?WT-EMP [54] [121] [54] 1-Water 98.3 99.5 99.9 98.6 99.9 99.9 2-Trees 91.2 91.7 95.7 93.5 97.8 97.4 3-Meadow 96.8 85.3 98.0 95.9 98.5 99.1 4-Bricks 88.4 99.2 99.7 98.8 99.7 99.8 5-B. Soil 97.0 98.4 99.8 99.4 99.8 99.9 6-Asphalt 96.3 98.6 98.7 98.4 99.0 99.5 7-Bitumen 96.0 98.1 97.1 98.2 98.3 98.5 8-Tiles 99.4 99.8 99.8 99.8 99.8 99.9 9-Shadows 99.9 99.2 99.9 99.9 99.9 100 OA 97.7 98.8 99.5 99.7 99.7 99.7 AA 95.6 96.6 98.8 98.1 99.2 99.3 κ0.96 – 0.99 0.98 0.99 0.99 Table 3.9: Classification results for the WT-EMP scheme on the Pavia City compared to SVM, EMP-ICA, and spec-EMP. The best accuracies are indicated in bold. The †indicates that the EMP is used as the only input to the SVM and ?that the stacked vector is built from the original data and the EMP. (a) (b) (c) Figure 3.13: WT–EMP on the Pavia City dataset: (a) true color representation, (b) spectral classification map with SVM, (c) WT-EMP classification map. 3.13 shows from left to right the true color representation, the spectral classification map with SVM and the WT-EMP classification map. 3.3. Scheme based on morphological profiles: WT–EMP 91 Indian Pines dataset The third scenario used in the experiments is the Indian Pines dataset, with the sixteen classes extracted from its reference map. This dataset has 220 spectral bands. A threshold of 0.01 for denoising the hyperspectral data obtained the best classification accuracy in this image. The spectral dimensionality is reduced 6 times by the 1D-DWT, creating stack vectors of 256 features (220 spectral bands + 36 features from the EMP). The number of samples used for training the SVM is the same as in [58, 108]. The training samples are randomly taken from the reference map, described in Section 2.9.3. Table 3.10 gives the classification accuracies obtained by the pixel-wise SVM classifier, by the WT–EMP scheme, and two MM-based schemes, KPCApand NW-NW [108]. The best results are indicated in bold. The MM-based schemes are based on EMAPs, creating the profile using the area and standard deviation attributes. The morphological profiles are created on the principal components extracted by kernel-PCA (KPCAp) and NWFE (NW-NW) feature extraction techniques, respectively. The first has 325 features that are used as input to the SVM classifier, and the NW-NW uses a feature extraction technique a second time to reduce the EMAP to 30 features. Two segmentation-based schemes, WSHED+MV and HSeg+MV, are also included for a broader comparison. The results of these two schemes are taken from [58] and not from [155], as the results presented by Fauvel were produced in [58] by the same number of training samples used in this experiment. (a) (b) (c) Figure 3.14: WT–EMP on the Indian Pines dataset: (a) true color representation, (b) spectral classification map with SVM, (c) WT-EMP classification map. 92 Chapter 3. Spectral-spatial classification schemes based on segmentation and MM SVM WSHED+MV Hseg+MV KPCApNW-NW WT–EMP [58] [58] [58] [108] [108] 1-Alfalfa 74.4 94.9 92.3 94.87 94.87 91.5 2-Corn-notill 78.2 94.2 90.5 71.41 90.24 83.7 3-Corn-mintill 69.6 78.1 83.0 92.73 98.85 85.7 4-Corn 91.9 88.6 95.7 96.20 97.28 94.6 5-Grass-pasture 92.2 95.1 94.4 93.96 95.52 92.3 6-Grass-trees 91.7 98.8 97.6 98.42 99.56 97.2 7-Grass-pasture-mowed 100 100 100 90.91 100 100 8-Hay-windrowed 97.7 99.5 99.5 98.63 99.54 99.1 9-Oats 100 100 100 100 100 100 10-Soybean-notill 82.0 96.3 92.1 87.39 86.27 84.4 11-Soybean-mintill 58.0 68.8 84.1 87.84 94.58 74.9 12-Soybean-clean 87.9 90.8 95.4 96.35 93.61 84.6 13-Wheat 98.8 99.4 98.2 99.38 99.38 99.2 14-Woods 93.0 99.7 98.6 95.90 92.36 87.7 15-Bldg-Grass-Trees-Drives 61.5 69.4 82.1 96.36 99.09 94.3 16-Stone-Steel-Towers 97.8 95.6 100 100 100 96.9 OA 78.2 86.6 90.8 90.20 94.2 86.6 AA 86.9 91.6 94.0 96.64 96.3 92.3 κ0.75 0.85 0.90 0.888 0.93 0.848 Table 3.10: Classification results for the WT–EMP scheme on Indian Pines and compared to SVM, WSHED+MV, HSeg+MV, KPCAp, NW-NW. The best accuracies are indicated in bold. 3.3. Scheme based on morphological profiles: WT–EMP 93 As can be seen from Table 3.10, the best OA (94.2%) and kvalue (0.93) are achieved by the NW-NW scheme, and the best AA (96.64%) for the KPCApapproach, both based on attribute profiles. Our proposed scheme gives an OA similar to the WSHED+MV scheme (86.6%), showing difficulties in classifying an image with low spatial resolution. However, the CS of the grass-pasture-mowed and oats classes are among the best, with a value of 100%. Even though, the EMP is considered an effective approach for combining spectral and spatial information [58, 108], the EMAP gives more accurate results in classifying the Indian Pines dataset. Figure 3.14, shows from left to right, the true color representation of the Indian Pines scene, the spectral thematic map with SVM, and the WT-EMP thematic map. Experimental results in presence of noise This section analyzes how each one of the stages of the proposed scheme works separately in the presence of noise. The two stages are WT–EMP‡, in which the Double-Density DWT is applied to the hyperspectral data, and the WT–EMP†stage, which only considers the EMP created from the wavelet coefficients as input to the SVM classifier. The experiments were carried out with the same configuration as described at the beginning of Section 3.3.4. The University of Pavia dataset was corrupted by additive white Gaussian noise (AWGN) [98] with a peak signal-to-noise ratio (PSNR) of 10, 16 and 20 dB. SVM WT–EMP‡WT-EMP†WT-EMP 1-Asphalt 56.5 90.2 89.4 97.6 2-Meadows 40.2 99.1 85.8 97.8 3-Gravel 15.9 99.6 83.9 98.4 4-Trees 78.5 80.5 91.9 96.4 5-Metal 91.9 99.8 99.5 99.6 6-Bare Soil 33.7 99.9 90.8 99.2 7-Bitumen 16.3 99.9 89.5 98.9 8-Bricks 49.2 98.2 85.8 97.4 9-Shadows 90.8 67.2 98.6 99.0 OA 45.9 96.0 88.0 97.3 AA 52.6 92.7 90.6 98.3 κ0.33 0.94 0.84 0.97 Table 3.11: Classification results for WT–EMP in presence of noise for a PSNR of 10 dB. University of Pavia dataset. 100 Chapter 4. Techniques and strategies for efficient GPU computing (a) (b) Figure 4.2: Hyperspectral data partitioning techniques: (a) spatial-domain partitioning, (b) spectral-domain partitioning. one byte, one value per bit, and the RGB values of a color image can be packed into three bytes. Thus, by packing data the number of global memory accesses is reduced when data are required in a kernel. In addition, if data from neighboring pixels are also required, the global memory accesses are even more reduced. This technique is used in the CA–Watershed based on Block–Asynchronous computation. 4.1.3 Spectral and spatial partitioning In hyperspectral imaging, two types of data partitioning techniques can be exploited: one in the spectral-domain and other in the spatial-domain [127]. These two data partitioned strategies are shown in Figure 4.2 where the [x,y]front plane in each figure corresponds to the first band of the image. In the spatial-domain partitioning, the pixel vectors are kept as a whole, as shown in Figure 4.2(a). In the spectral-domain partitioning, the image is subdivided into slices comprising contiguous spectral bands, as shown in Figure 4.2(b). The data partitioning approach strongly depends on the processing techniques applied to analyze the hyperspectral image. For example, if all the spectral features of a pixel are required at once to make a computation, as in the case of unmixing processing, the spatial-domain partitioning is desired, as each pixel vector can be assigned to a different block of threads. However, if the computation can be individually done in each band, as in the case of morphological operations, 4.1. Introduction 101 the spectral-domain partitioning can take advantage of the fine-grained level of parallelism of CUDA. In this work, different hyperspectral data partitioning strategies and thread block arrangements are studied in order to effectively exploit the memory and computing capabilities of the GPU architecture. 4.1.4 Challenges of GPU computing In order to maximize the performance on the GPU, different aspects must be taken into consideration. This section highlights some key GPU performance issues that can be grouped in three main rules for GPGPU programming: 1) minimizing data transfer between CPU and GPU, 2) giving the GPU enough work to do, and 3) focusing on data reuse within the GPU [53]. Most of the important issues are related to an efficient use of the memory hierarchy. Data movement from the CPU to the GPU global memory is through the PCIe bus, which is very slow (8GB/s) in comparison to the peak bandwidth of the global memory (250 GB/s) and the shared memory (2.5 TB/s). Therefore, data should be copied to the GPU once and be reused as much as possible. As memory operations are executed per warp, it is important to known how data are going to be accessed, which mainly depends on how we implement the algorithm. By aligning accesses to consecutive memory locations in global memory, we ensure coalesced accesses (see CUDA parallel programming model in Section 2.8.2), minimizing load/store operations into the fewest possible memory transactions. Data in constant memory can be broadcast if the same value is accessed by all the threads within the warp. And texture memory has special features, such as filtering, and a dedicated cache memory which can improve performance when threads access values in some regular spatial neighborhood, for example in a 3 ×3 window. Therefore, data may need to be rearranged or distributed in a different way in order to obtain maximum performance of the GPU memory. Data packing, as described in Section 4.1.2, can be used to minimize data transfers, and different data partitioning strategies, see Section 4.1.3, as well as thread block arrangements can improve the performance of the GPU. The size of a block of threads should be a multiple of the warp size (32). In most of the cases, it is often necessary to perform different kernel launch configurations (threads per block and blocks per grid) to find the best one that maximizes the device utilization [118]. 102 Chapter 4. Techniques and strategies for efficient GPU computing The limit in the hardware resources required to execute a kernel can be achieved by the size of the block, as described in Section 2.9.1. The occupancy, which can be expressed as the number of concurrent threads per SM can be used to measure the performance on the GPU. The amount of occupancy required to reach the maximum performance depends on the code. If the code is limited by memory accesses, the highest number of concurrent blocks per SM is desired in order to hide the latency of those accesses. Instructions are also executed per warp, so the 32 threads within the warp take the same path in the code. In conditional control flows, such as if/else statements, if at least one thread within the warp takes a different path, we have a warp divergence in the code, which potentially affect performance [118]. Hence, care must be taken to write a coherent control flow code. In some cases, sequential code must be rewritten to expose sufficient parallelism to the GPU, and arithmetic instructions reordered and split for balancing the workload among memory accesses and computation. Regarding data reuse within the GPU, if data accesses have sufficient locality, redundant accesses to global memory must be minimized by exploiting the shared and the cache memories. The main use for the shared memory is to reuse data within a block and share data among the threads of the same block. This memory is managed by the programmer, and the effective use of this memory can lead to speedups of 7×compared to global memory implementations [40]. This is a significant challenge when programming for the GPU as data in this on-chip memory are only shared for the threads of the same block. Sharing data among all the threads must be through the global memory, coming into play the aforementioned key GPU performance issues related to an efficient use of the memory hierarchy. Atomic memory operations may be necessary to avoid race conditions in applications where multiple threads need to access the same memory space simultaneously for reading and writing data, and there is no other way of synchronizing. CUDA provides atomic operations for ensuring that all concurrent updates to the same memory location can be performed atomically, so that all threads can see those updates. However, the atomic operations can degrade the performance if many threads try to carry out an atomic operation on a small number of memory locations. A extensive range of skills are needed to be an effective parallel programmer [90]. These skills can be grouped into computer architecture, programming models and compilers, algorithm techniques and domain knowledge. The CUDA C programming guide [119] and the 4.2. Block–Asynchronous strategy 103 CUDA C best practices guide [118] have a comprehensive manual and numerous tips and key issues for increasing the computational throughput of NVIDIA GPUs. 4.2 Block–Asynchronous strategy In this section we present the Block–Asynchronous (BA) strategy proposed in this thesis for efficient computation of problems which require iterative computation, such as the Cellular Automata (CA). This method reduces the number of points of global synchronization allowing efficient exploitation of the memory hierarchy of the GPU. By BA computation, we mean updating a group of values an unbounded number of times without a global synchronization. Thus, each region is updated asynchronously with respect to other regions. For example, a cellular automata can be partitioned into different groups of cells which can be updated independently and locally. This strategy perfectly matches the tiling/grids parallel pattern described in the preceding section. This model is shown in Fig. 4.3 for a 6×6 cellular automaton. Each block (tile) of 3 ×3 cells is updated an unbounded number of times (intra-block updates) but the values outside the block, corresponding to the apron, are kept constant; i.e., equal to their values at the beginning of the stage. The entire grid is updated after a global synchronization, so data are read at the block boundaries (inter-block updates), which allows the propagation of data across the blocks. The BA strategy is also adequate for multicore architectures. This strategy can easily be adapted to solve different algorithms, such as the asynchronous cellular automaton to compute the watershed transform on GPU [139, 135, 138]. Advanced MM operations can also benefit from the BA approach, for example opening and closing by reconstruction and greyscale attribute filtering, as it will be shown in the next sections. The remainder of this section is organized as follows. Section 4.2.1 present the CA– Watershed algorithm based on Block–Asynchronous computation and its extension to 3D for volumes processing. In Section 4.2.2 we present a Block–Asynchronous algorithm for computing opening and closing by reconstruction (BAR) on GPU [134]. In Section 4.2.3 we describe a new proposal for greyscale area opening and closing on GPU which can be extended to other attributes. Finally, the results are discussed in Section 4.2.4. 104 Chapter 4. Techniques and strategies for efficient GPU computing Figure 4.3: Example of Block–Asynchronous strategy, mapping on 6×6 CA (4-connectivity). Intra-block updating (left) and inter-block updating (right). 4.2.1 CA–Watershed based on Block–Asynchronous computation on GPU The Watershed Transform based on cellular automata (CA–Watershed) was described in Section 3.2.2. The asynchronous behaviour of this automaton, introduced by Galilée et al. [64], is subject to its implementation. In the following, we describe first the synchronous implementation that can be executed in CPU using OpenMP, as well as on GPU. The idea is giving the foundations for understanding the asynchronous behaviour of the automaton. Then, the GPU Block–Asynchronous computation of the same algorithm is described. The algorithm is non-deterministic and may introduce artifacts into the border of the segmented regions. Therefore, we propose an artifacts-free block–asynchronous GPU implementation that produces the correct results by correcting the data propagation speed among the blocks using wavefront techniques [112]. These implementations follow the parallel tiling/grids pattern in which the grid of cells of the automaton is partitioned into regular regions that are assigned to blocks of threads of the GPU. Finally, due the nature of the CA, the implementation can be easily extended to three dimensions in order to process a 3D volume. Block–Synchronous computation on GPU The CA–Watershed synchronous implementation has two kernels: one for initializing and another for updating the automaton. The pseudo-code presented in Figure 4.5 shows the 4.2. Block–Asynchronous strategy 105 Compute and non-minimum MP NM INIT look at look at Extend plateaus Hill-Climbing Figure 4.4: Three-state cellular automaton implementing Hill-Climbing algorithm [64]. block-synchronous computation on GPU. The kernels executed on GPU are placed between <> symbols. The pseudocode also includes the Tex, GM and SM acronyms to indicate kernels executed with data on texture, global and shared memory, respectively. Figure 4.4 in this page recalls the 3–state cellular automaton described in Section 3.2.2 to compute the watershed transform. The NF(u)and the set of neighbors with the same grey value, N=(u), are computed for each pixel in the first kernel (line 2 in the pseudocode). In the second kernel (line 4), the updates flood each region with a representative label in an iterative process executed by the CPU with a global synchronization at each step (line 5). These kernels are configured to work in rectangular thread blocks with a thread operating on one pixel. With the first kernel, the automaton is initialized according to l(u) = min(l(v)) and f(u) = fNF(u), corresponding to Eq. (3.1) and Eq. (3.2), respectively (see page 68). The grey values are read from texture memory, as this read-only memory speeds up the ac- Input: one band image Ion the global memory of the GPU Output: segmentation map 1: copy input data Ifrom global to texture memory 2: <initialize cellular automaton> . (Tex and GM) 3: while cellular automaton is not stable do 4: <synchronous updating of the automaton> . (GM) 5: global synchronization among blocks 6: end while Tex: texture memory, GM: global memory. Figure 4.5: Pseudocode for CA–Watershed synchronous implementation on GPU. 106 Chapter 4. Techniques and strategies for efficient GPU computing (a) (b) Figure 4.6: Example of data packing. (a) Data of the automaton packed into 64 bits, (b) structure of the variable N=(u)using 4-connectivity packed in 1 byte. cesses to data when they present high spatial locality. Once all data have been initialized, they are packed into 64 bits before being transferred to global memory, with the objective of increasing the data locality and reducing number of global memory accesses. The minimum amount of data required are 4 bytes for l(u), 1 byte for h(u), 1 byte for N=(u)and 2 bytes for NF(u). Fig. 4.6(a) shows an example of data packing for one pixel. The set N=(u)is compressed in 1 byte using 1 bit for each neighbor (L, R, U, D in Fig. 4.6(b)) where “1” means a neighbor with the same grey level and “0” means a neighbor with a different one. The four least significant bits are ignored but they may be used to connect up to 8 neighbors. It is not necessary to store the state of each pixel as it can be deduced from NF(u)(see Figure 4.4). If the lower slope is empty, the state is MP; otherwise, the state is NM. At the end of the initialization stage the state of each pixel, a cell of the automaton, has switched to NM or MP. In order to update all the pixels synchronously, we use two 64-bit buffers, one input buffer for reading data and one output buffer for writing the results. The updating stage has been implemented through a loop executed by the CPU (lines 3–6 in the pseudocode of the Figure 4.5), which calls a CUDA kernel at each step (line 4). There is one global synchronization per step, as indicated in line 5, which makes the computation wait until the GPU finishes its current work [119]. The input and output buffers are swapped before the next iteration. Only one flag needs to be moved to the CPU at each inter-block iteration 4.2. Block–Asynchronous strategy 107 Input: cellular automaton packed in 64 bits Output: cellular automaton updated 1: load and unpack data from global memory to registers .input buffer 2: updates image and labels according to Eqs. (3.3)–(3.5) 3: pack results in 64 bits from registers and store data in global memory .output buffer Figure 4.7: Pseudocode for CA–Watershed synchronous CUDA kernel executed in global memory. indicating whether the automaton must be further processed. The kernel for the synchronous updating of the automaton is shown in Figure 4.7. In each call to the kernel, data are read from the input buffer in global memory and are unpacked and automatically stored in registers (line 1). The pixels are updated once (line 2) as described by Eq. (3.3) and Eq. (3.4) if their state is MP, and by Eq. (3.5) if their state is NM. These equations are described in Section 3.2.2. Finally, the resulting data are packed and stored in the output buffer. The update ends when all regions have been flooded. Block–Asynchronous computation on GPU The CA–Watershed based on block-asynchronous computation has the advantage of reusing information within a block, unlike the synchronous one, efficiently exploiting the shared and cache memories of the GPU. The storage requirements are the same as for the blocksynchronous implementation. The updating stage has been adapted to perform in shared memory as many updates inside a region as possible (intra-block updates), before performing a synchronization among thread blocks (inter-block updates). Each region is synchronously updated (i.e. all cells within a region are updated at each iteration), while the regions themselves are asynchronously updated (an update of the entire grid is performed at certain selected steps). Hence, this is a hybrid iterative process that includes block-asynchronous and blocksynchronous updates. The pseudocode in Figure 4.8 and Figure 4.9 show the inter-block updates (lines 3–6) and the intra-block updates (lines 2–5), respectively. The kernels is placed between the <and >symbols. The inter-block updating loop is executed on the CPU and calls the asynchronous updating kernel that is executed on the GPU. In this kernel, for each block, once data are loaded in shared memory, the pixels are modified according to Eqs. (3.3) – (3.5) in an iterative intrablock process within each region of the image. Threads within a block are synchronized 108 Chapter 4. Techniques and strategies for efficient GPU computing Input: one band image Ion the global memory of the GPU Output: segmentation map 1: copy input data Ifrom global to texture memory 2: <initialize cellular automaton> . (Tex and GM) 3: while cellular automaton is not stable do .inter–block updating 4: <asynchronous updating of the automaton> . (SM) 5: synchronization among blocks .global synchronization 6: end while Tex: texture memory, GM: global memory, SM: shared memory Figure 4.8: Pseudocode for Asynchronous CA–Watershed implementation on GPU. Input: cellular automaton packed in 64 bits Output: cellular automaton updated 1: load and unpack data from global memory to shared memory .input buffer 2: while cellular automaton is not stable within a block do .intra–block updating 3: updates image and labels according to Eqs. (3.3)–(3.5) .(SM) 4: synchronize threads within the block .local synchronization 5: end while 6: pack results in 64 bits from shared memory and store data in global memory .output buffer Figure 4.9: Pseudocode for Asynchronous CA–Watershed CUDA kernel executed in shared memory. locally at each step of the intra-block process, so data updated within a block can be reused from the shared memory, which is much faster than the global memory space (see Section 4.1.4). The intra-block updating ends when no new modifications are made with the available data within the region. Then the data in shared memory are packed and stored in global memory. The updating stage ends when all regions have been flooded. In order to update the pixels at the edge of the block in this block-asynchronous implementation, the shared memory allocated for each region must be extended with a border of size one. This extra shared memory space corresponds to the apron illustrated in Figure 4.1(b) that is required in tiling/grid parallel patterns when the local computation also involves neighboring data. Thus, the apron of one region overlaps the adjacent regions. Threads on the edge of the block have to do extra work loading the data of the border. 4.2. Block–Asynchronous strategy 109 The block-asynchronous algorithm avoids points of global synchronization that are costly in execution time, and efficiently exploits the shared memory, which has lower access times than the global memory. Artifacts-free block–asynchronous computation on GPU The block–asynchronous CA–Watershed implementation obtains a correct segmentation according to the watershed segmentation definition. Thus, when non-minimum plateaus exist in the image, the algorithm gives a correct segmentation; nevertheless, the watershed lines may not match the geodesic distance properly. When the data propagation speed is similar for all the cells, the algorithm may give a good approximation of the watershed lines. However, when computed by regions, it presents the problem of data propagation at the region boundaries, which causes artifacts as shown in the example in Figure 4.10. Initially data are propagated within the region during the intra-block updating, and later this is performed at the region boundaries during the inter-block updating. The different speed of data propagation in the intra-block and inter-block updates result in improperly placed watershed lines. This situation is visually observed as small irregularities in the watershed lines, as illustrated in Figure 4.10. As a consequence of the asynchronous computation by blocks some undesirable artifacts arises and the quality of the segmentation is slightly affected. The asynchronous behavior of the CA–Watershed implementation is fixed by correcting the data propagation speed among the blocks. The artifact–free block-asynchronous algorithm is based on the application of a technique known as wavefront [112] increasing the quality of the watershed lines obtained. The wavefront starts from the lower border of a plateau (see Section 2.3.2) and iteratively Figure 4.10: An image of size 128×128 pixels (left), the correct watershed line obtained using 4-connectivity (middle), and the artifacts produced by the asynchronous computation using blocks of cells of size 32×32 (right). 116 Chapter 4. Techniques and strategies for efficient GPU computing (a) (b) Figure 4.14: Example of a greyscale attribute opening. (a) Greyscale image with five connected components (the grey level is indicated by the subscript, and the i-th component within a level is indicated by the superscript) and the corresponding max-tree, (b) result of the attribute opening on the greyscale image and the pruning of the max-tree. can be constructed by merging the different sub-trees using the unique labels assigned by the union-find algorithm. Direct implementations without building the max-tree are also possible based on hierarchical FIFO and priority queues [170, 43, 97]. In [43] the filtering and the flooding for labelling the components are combined by passing the max-tree construction. It is considered a direct approach only for area filtering. The work presented in [97] is based on [43] but it can 4.2. Block–Asynchronous strategy 117 Input: greyscale image I Output: attribute opening 1: copy input data Ifrom CPU to the global memory of the GPU 2: <labelling the connected components of Iat each grey level> . BA labelling (SM) 3: for each grey level h∈histogram(I)do 4: <merge the regions of the connected components at level h> . (GM) 5: <calculate the attribute for each region> . (GM) 6: <perform attribute filtering> . (GM) 7: end for GM states for computation in global memory and SM in shared memory. Figure 4.15: Pseudocode for the greyscale attribute opening on GPU. be applied to attributes other than area. One drawback of these approaches that do not create the max-tree is that the pruning rule must be known a priori. The grey value assigned to the final image after removing a component cannot be retrieved by traversing the tree as it is not constructed. In this thesis we propose a GPU implementation for greyscale attribute openings and closings without using queues, through threshold decomposition and merging of connected components. This proposal is a member of the group of direct implementations that simulate the max-tree [170, 43, 97]. The pseudocode in Figure 4.15 shows the workflow for the greyscale attribute opening on GPU proposed in this thesis. The kernels executed on GPU are placed between <> symbols. The area closing can be computed with the same algorithm but using the complement of the image. The connected component labelling (line 2) is performed by using the block–asynchronous strategy proposed in this thesis. All the connected components of the image are labelled at once. First, each pixel is given a unique label that identifies its position within the image in a row-major order. Then, an iterative process propagates the minimum label between all the connected neighbors. This is similar to the block–asynchronous computation described on page 108. The updating stage is adapted to perform in shared memory as many intra-block updates as possible, before performing the inter-block synchronization. Second, the image is processed through threshold decomposition (lines 3–7). At each grey level h, from the highest to the lower intensity, the connected components at grey level h, which have been already labelled in the last iteration, are merged if they have a common 118 Chapter 4. Techniques and strategies for efficient GPU computing border (line 4). This iterative process simulates the max-tree from the leaves to the root and only keep in memory the nodes at the current level h. Once the components have been merged, the attribute for each one can be computed (line 5) in a new kernel. Finally, based on the value of the attribute, a component is filtered at the current level h (line 6) if the criterion for that component is false. We have used the filtering max rule that prunes the branches from the leaves up to the first node that needs to be preserved [144]. This proposal could be used to compute other attributes than the area of a region, simply by modifying the kernel at line 5. 4.2.4 Results In this section we present the performance results for the Block–Asynchronous strategy applied to the asynchronous cellular automaton to compute the watershed transform, the opening and closing by reconstruction and the area attribute filtering. Experimental setup The different implementations based on the Block–Asynchronous (BA) strategy have been evaluated on the Intel quad-core i7-860 microprocessor. The GPUs used in the experiments are the GTX 580 and the GTX TITAN based on the Fermi and Kepler architectures, respectively. The CUDA code has been compiled under Linux using the nvcc compiler with the CUDA toolkit 4.2 (GTX 580) and 5.5 (GTX TITAN). The hardware, the compute capability of the graphic cards and the images used in these tests are described in Section 2.9.1. The reference codes for comparison are optimized OpenMP parallel implementations for the CA–Watershed, the Fast Hybrid Reconstruction (HRA) algorithm [169] for the opening and closing by reconstruction, and an efficient implementation based on the max-tree (mintree) [42] for attribute filtering. These algorithms will be described in the subsection where the comparisons are made. The performance results analyzed are expressed in terms of execution times and speedups. The execution times were obtained as the average of 20 executions. In all the tests, a connectivity of four pixels is used. The datasets used in the experiments are two images (Lena and CT Scan Head) and the BrainWeb volume. The datasets are described in Section 2.9.3. The two images used are representative cases of processing small (Lena) and large (CT Scan Head and BrainWeb) plateaus, respectively. The processing of large plateaus (regions of uniform grey values) makes it necessary to propagate the labels through large regions of the image. This requires 4.2. Block–Asynchronous strategy 119 Size 512×512 1024×1024 2048×2048 Transfer time 0.0022s 0.0082s 0.0321s Size 45×54×45 90×108×90 181×217×181 Transfer time 0.0013s 0.0073s 0.0581s Table 4.1: CPU–GPU data transfer times for 2D and 3D images at different sizes. more computation time than processing small plateaus. Thus the selected images represent two very different cases regarding to computational cost of data propagation among blocks. The CPU–GPU data transfers are carried out at the beginning and are shown in seconds in Table 4.1. This time is the same independently of the application where the BA approach is used. CA–Watershed based on Block–Asynchronous computation on GPU In this section we present the results for the following GPU implementations of the watershed transform based on cellular automata: block-synchronous, block-asynchronous and artifactsfree block-asynchronous implementations1. First, we have checked the correctness of the asynchronous CA–Watershed implementation by comparing the number of segmented regions obtained by the GPU algorithms to the number of regions obtained by a sequential watershed algorithm over the images. Table 4.2 shows the number of regions created by the watershed transform, which is the same for the CPU and GPU implementations. The difference between the number of regions at different resolutions is due to the process of scaling the image. The CT Scan Head and the BrainWeb datasets present large plateaus and therefore the number of regions is lower (but larger in size) than for the Lena image. In the following, the performance is analyzed in terms of occupancy, execution times and speedup. – Analysis of the GPU Parameters: The initialization and updating kernels used on the GPU implementation, see the pseudocode presented in Figure 4.7 and Figure 4.9, have been 1Part of these results have been published in P. Quesada-Barriuso, D. B. Heras, and F. Argüello, “Efficient 2D and 3D watershed on graphics processing unit: block-asynchronous approaches based on cellular automata,” Computers & Electrical Engineering, vol. 39, no. 8, pp. 2638–2655, 2013. 120 Chapter 4. Techniques and strategies for efficient GPU computing Size 512×512 1024×1024 2048×2048 Lena 24958 25139 28521 CT Scan Head 6221 7300 13381 Size 45×54×45 90×108×90 181×217×181 BrainWeb 1115 5669 14348 Table 4.2: Number of regions generated by the watershed transform. analyzed according to the resources available on the GTX 580 (Fermi architecture). This GPU has 16 SMs with the following limits per SM: 1536 threads, 8 active blocks, 32768 registers of 32 bits and 64 KB of on-chip memory that can be configured as a shared memory of 16 KB and 48 KB for the L1 cache or vice versa (see Table 2.3 for a full description of the GPU). Table 4.3 shows the maximum number of active blocks based on the block size and the on-chip memory configuration. For the ca_watershed_asynchronous kernel, the reason for selecting the L1 configuration with the remaining 16 KB being for the shared memory is that the maximum number of active blocks is given by the limit of 1536 threads per SM. As described in Section 4.2, 8 bytes per pixel are required for data packing. For a block with 16×16 threads, the shared memory required would be 2048 bytes but as the shared memory has been extended with an apron of size one, each block needs 18×18×8 bytes of this onchip memory, i.e. 2.5 KB per block. Therefore, considering a maximum of 6 blocks (1536 threads) per SM, a total of 15 KB of shared memory per SM are used. For the case of the artifacts-free asynchronous proposal (ca_watershed_asynchronous- free in Table 4.3), the number of simultaneously active blocks per SM is reduced to four. In this case the limiting factor is the number of 32768 registers available per SM. The proposal requires 26 registers per thread, which gives a total of 16×16×26 =6656 registers per block, Kernel 16×16 32×16 32×32 8×8×4 ca_watershed_asynchronous L1 6 3 1 3 ca_watershed_asynchronous-free L1 4 2 1 2 ca_watershed_asynchronous Sh 6 3 1 6 ca_watershed_asynchronous-free Sh 6 3 1 6 Table 4.3: Number of active blocks per SM for the different kernels based on the block size and the shared memory requirements. L1 indicates that 48 KB are used for the L1 memory and 16 KB for the shared memory. Sh states for the opposite configuration. Analysis for CA–Watershed implementations. 4.2. Block–Asynchronous strategy 121 so there are enough registers in each SM for only 4 blocks. Regarding the shared memory use in the artifacts-free kernel, each pixel requires 4 extra bytes to manage the geodesic distance properly, (see Algorithm 2), so 18 ×18×(8+4) = 3888 bytes per block are required. By using the Sh configuration (48 KB for shared memory), the limit in the number of active blocks is given by the block size. – Performance analysis: The OpenMP implementation is based on the block-synchro- nous approach (see pseudocode on Figure 4.5) and it uses 4 threads scheduling the work statically among the threads by a loop construct distributing the iterations into 4 chunks of the same size (one per thread), in order to evenly distribute the workload among the threads and to achieve a high locality in the data accesses. The need to access data outside the region assigned to each thread is not a problem in the OpenMP implementation, as all the threads access the same memory space. The algorithm includes an implicit synchronization barrier at each step of the updating. Table 4.4 gives the performance results obtained in the GTX 580. The execution times for the GPU proposals in this table include the CPU–GPU data transfer time. In all the tests, the CA–Watershed based on block-asynchronous computation obtains high speedups for all the image sizes. As shown in Table 4.4, the speedups also scale well with the size of the image; i.e. from 9.0×to 13.6×for the block–synchronous implementation for the Lena image up to 11.7×to 22.2×with the block-asynchronous proposal. When the image size increases, so does the amount of computational work, the hundreds of available threads are better exploited. The 3D GPU proposals obtain speedups for all the volume sizes, and the speedup values increase with the volume size as the computational load also increases. The performance results for the block-asynchronous proposals are always better than for the synchronous implementation; approximately twice as good. The performance results for both, the block-asynchronous and the artifacts-free block-asynchronous approaches are very similar. We focus now the test on the 2D images at a resolution of 2048×2048 pixels as the behavior of processing large plateaus is better appreciated in this case. When the image presents large plateaus the computational cost of the watershed transform increases as the labels must be propagated through large regions of the image. Comparing the execution times of the block-synchronous and block-asynchronous approaches (see Table 4.4), a speedup of 1.6x is obtained for the image of Lena while the speedup increases up to 4.3x for the CT Scan image. The improvement of the block-asynchronous proposal versus the synchronous imple- 122 Chapter 4. Techniques and strategies for efficient GPU computing Lena (2D) 512×512 1024×1024 2048×2048 OpenMP (4 threads) 0.0351s 0.1990s 1.2452s GPU Synchronous 0.0039s (9.0×) 0.0188s (10.6×) 0.0916s (13.6×) GPU Asynchronous 0.0030s (11.7×)0.0131s (15.2×)0.0562s (22.2×) GPU Artifacts-Free Async. 0.0034s (10.3×) 0.0158s (12.6×) 0.0736s (16.9×) CT Scan (2D) 512×512 1024×1024 2048×2048 OpenMP (4 threads) 0.4941s 2.8793s 15.0919s GPU Synchronous 0.0305s (16.2×) 0.1436s (20.1×) 0.6992s (21.6×) GPU Asynchronous 0.0093s (53.1×)0.0381s (75.6×)0.1628s (92.7×) GPU Artifacts-Free Async. 0.0126s (39.2×) 0.0522s (55.2×) 0.2353s (64.1×) BrainWeb (3D) 45×54×45 90×108×90 181×217×181 OpenMP (4 threads) 0.1337s 2.2611s 37.7378s GPU Synchronous 0.0084s (15.9×) 0.0820s (27.6×) 1.1227s (33.6×) GPU Asynchronous 0.0044s (30.4×)0.0451s (50.2×)0.5907s (63.9×) GPU Artifacts-Free Async. 0.0050s (26.5×) 0.0540s (41.8×) 0.7304s (51.7×) Table 4.4: Performance results including data transfer times (speedup in brackets). Best results in bold. mentation is better for the second image; although, as shown in Table 4.4, processing large plateaus takes more time: 0.1628s for the CT Scan image while the Lena image only requires 0.0562s. This is because for the block-asynchronous proposal the intra-block updating allows the labels to propagate faster among regions, especially in images with large plateaus. If a region is entirely within a plateau, the labels have to be propagated from side to side of that region. In this situation, only one inter-block update and wintra-block updates are needed, where wis the width of the region. The synchronous implementation would need winterblock updates, with the consequent penalty for transferring data from and to global memory at each step, with each one of those steps corresponding with a global synchronization. The block-asynchronous approach reduces the number of synchronizations among thread blocks and increases data reuse thanks to the inter- and intra-block updating scheme. The decrease in the number of synchronizations for the block-asynchronous proposals is illustrated in Table 4.5, where the number of inter-block and intra-block updates are summarized for the synchronous and the block-asynchronous implementations and the test images. Only the values for the block-asynchronous implementation are shown, as the numbers are the same for the artifacts-free proposal. For the block-synchronous implementations (on CPU and on GPU) only inter-block updates take place in the sense that after each update of all the 4.2. Block–Asynchronous strategy 123 Lena inter-block intra-block (min.) intra-block (max.) intra-block (avg.) GPU Synchronous 114 — — — GPU Asynchronous 16 22 195 108.5 CT Scan Head inter-block intra-block (min.) intra-block (max.) intra-block (avg.) GPU Synchronous 1156 — — — GPU Asynchronous 76 88 1248 668 Table 4.5: Number of updates per pixel for the block-synchronous and block-asynchronous implementations for the 2048×2048 images. pixels of the image one global synchronization operation is required. Observing, for example, the values for the CT Scan image in the table, the number of inter-block updates (i.e. the number of global synchronizations required) is 1156 for the synchronous implementation. For the block-asynchronous cases the number of inter-block updates decreases to 76 and the total number of asynchronous intra-block updates per block summing up all the iterations ranges from 88 to 1248, depending on the block, with 668 being the average value over all the blocks. Hence, the number of updates per pixel is 1156, with the same number of corresponding global synchronizations for the synchronous implementation, and an average of 668 local synchronizations with only 76 global synchronizations for the block-asynchronous algorithms. In the case of the Lena image, a similar decrease is observed. – Comparison to other works: Proposals of different watershed algorithms on the GPU using shaders [88] and CUDA [171, 91, 81] have been presented in the last few years. In [171] a new algorithm is presented based on the introduction of a chromatic function for establishing the order in which the voxels are processed. The experiments are carried out over volume data sets, obtaining maximum speedups of 7×on a GTX295 when compared to the sequential proposal of the algorithm, even when large volume data sets of up to 600 ×600 ×600 are considered. In our case, the largest volume considered was 181×217×181, 9 times smaller, achieving speedup values of 63.9×on a GTX 580. Taking into account that the speedups of our block-asynchronous proposals increase with the volume size, as shown in Table 4.4, and that our experiments have proved that the block-asynchronous proposals scale by a factor of 2×to 4×between the GTX 295 and the GTX 580 GPUs [139], we can conclude that our proposals outperform the results in [171]. 124 Chapter 4. Techniques and strategies for efficient GPU computing The algorithm presented in [91] is inspired by the drop of water paradigm and performs a component labelling and a path compression approach [74]. Its results are compared to a GPU synchronous algorithm for the watershed based on a CA [88], outperforming it. Given that the experiments in [91] are performed on an older GPU than in our case, we have executed them on our GTX 580, obtaining similar speedup results to the obtained with our blockasynchronous proposal of the CA-watershed described in this work. Opening and closing by reconstruction on GPU The opening and closing by reconstruction GPU implementation based on the block-asyn- chronous strategy, see Section 4.2.2, has been evaluated on the GTX TITAN (Kepler architecture)2. The block size is configured with 32 ×8 threads. Each block requires only 340 bytes of shared memory, so the on-chip memory is configured with 48 KB for L1 cache. The reference codes for comparison are the Fast Hybrid Reconstruction (HRA) algorithm in CPU [169], and the GPU Sequential Reconstruction (SR_GPU)3algorithm proposed in [87]. The performance results analyzed are expressed in terms of execution times and speedups. The speedups are calculated with respect to the HRA (CPU) and the SR_GPU algorithms. The execution times were obtained as the average of 20 executions. The datasets used in the experiments are the Lena and the CT Scan Head images. 2Part of these results have been published in P. Quesada-Barriuso, F. Argüello, D. B. Heras, and J. A. Benediktsson, “Wavelet-based classification of hyperspectral images using extended morphological profiles on graphics processing units,” Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, vol. PP, no. 99, pp. 1–9, 2015 (published online, print edition pending). 3My thanks to Prof. Karas for sharing his implementation. Lena 512×512 1024×1024 2048×2048 HRA (CPU) 0.1100s (1.0×) 0.1244s (1.0×) 0.1723s (1.0×) SR_GPU (GPU) 0.0168s (6.5×) 0.0610s (2.0×) 0.2109s (0.8×) BAR (GPU) 0.0027s (40.1×)0.0112s (11.1×)0.0485s (3.5×) CT Scan Head 512×512 1024×1024 2048×2048 HRA (CPU) 0.1086s (1.0×) 0.1243s (1.0×) 0.1541s (1.0×) SR_GPU (GPU) 0.0194s (5.6×) 0.0679s (1.8×) 0.2097s (0.7×) BAR (GPU) 0.0026s (41.7×)0.0113s (11.0×)0.0388s (3.9×) Table 4.6: Block–asynchronous morphological reconstruction results including data transfer times (speedup in brackets). Best results in bold. 4.2. Block–Asynchronous strategy 125 For this experiment we have created a marker image Jfor each dataset Ias J(p) = max{I(p)−h,0}, that is known as the hmax transform. A value h=10 was used to create the initial marker. The results for the three algorithms are shown in Table 4.6. Both, SR_GPU and BAR algorithms include the data transfer from CPU to GPU memory. It can be observed that the GPU proposals outperform the HRA algorithm. However, the SR_GPU algorithm cannot beat the HRA algorithm with the images of 2048×2048 pixels. It can be observed that the speedup decreases by increasing the size of the images. The reason is that the HRA algorithm is optimized on CPU to process only those pixels that need to be reconstructed, while the SR_- GPU and the BAR algorithms are based on a raster scanning of all the pixels of the image. The proposed algorithm offers a better performance in all cases. By performing multiple scans in shared memory at each iteration, this morphological reconstruction on GPU efficiently exploits the shared memory through the block-asynchronous updating process. The best speedup is 41.7×obtained on the CT Scan Head (512×512 pixels). Greyscale attribute filtering on GPU This section presents the results obtained by the GPU greyscale attribute opening (closing) algorithm, described in Section 4.2.3. We have evaluated the performance on the GTX TITAN (Kepler architecture). The available resources for this GPU are described in Section 2.9.1. The reference code for comparison is an efficient implementation based on the max-tree (min-tree) presented in [42] for classification of hyperspectral images by using extended attribute profiles (EAPs). The University of Pavia dataset, described in Section 2.9.3, is used in the test. The first principal component extracted by PCA from this dataset is used to create an EAPabased on area and an EAPdbased in the diagonal of the bounding box. The thresholds (λ) for creating each profile are λ∈{36,49,169,361,625,1369}for the area and λ∈{10,25,50,100,150,250}for the diagonal. These values were extracted from [41]. Max-tree / min-tree Opening Closing TOTAL (CPU) [42] GPU GPU GPU EAPa0.9450 0.0825 0.1200 0.2025 (4.6×) EAPd1.0331 0.0941 0.1335 0.2276 (4.5×) Table 4.7: Performance results for the greyscale attribute filtering (opening and closing) including data transfer times for the first PCA of the hyperspectral image of the University of Pavia.