Identificación de tumores cerebrales mediante el estudio de la forma y la composición de los tejidos usando imágenes hiperespectrales

UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA ESCUELA DE INGENIERÍA DE TELECOMUNICACIÓN Y ELECTRÓNICA PROYECTO FIN DE CARRERA Identificación de tumores cerebrales mediante el estudio de la forma y la composición de los tejidos usando imágenes hiperespectrales TITULACIÓN: Ingeniero de Telecomunicación TUTOR/ES: Gustavo Iv á n Marrero Callic ó Inmaculada Garc í a D ó pido AUTOR: Miguel Ángel Tejedor Hernández FECHA: Septiembre 2015 UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA ESCUELA DE INGENIERÍA DE TELECOMUNICACIÓN Y ELECTRÓNICA PROYECTO FIN DE CARRERA Identificación de tumores cerebrales mediante el estudio de la forma y la composición de los tejidos usando imágenes hiperespectrales Presidente: Tutores: Secretario: Vocal: Autor: NOTA: TITULACIÓN: Ingeniero de Telecomunicación TUTOR/ES: Gustavo Iv á n Marrero Callic ó Inmaculada Garc í a D ó pido AUTOR: Miguel Ángel Tejedor Hernández FECHA: Septiembre 2015 RESUMEN La presente memoria resume el trabajo de investigación realizado por Miguel Ángel Tejedor Hernández con motivo de su proyecto fin de carrera (PFC). En concreto, el presente trabajo describe una comparativa de diferentes cadenas de procesamiento de imágenes hiperespectrales de tumores cerebrales humanos. El documento sigue la estructura clásica de un trabajo de investigación en dicho campo, presentando en primer lugar las motivaciones y objetivos que han motivado la comparativa de diferentes técnicas de análisis hiperespectral respondiendo a una necesidad claramente existente en este campo de estudio, pues actualmente no existen mecanismos fiables para delimitar el tejido tumoral con el fin de extraerlo en su totalidad y de forma precisa durante una operación de neurocirugía. A continuación se realiza un estudio en profundidad del estado del arte en dicho campo, desde el concepto de píxel hiperespectral hasta los algoritmos existentes que fundamentan la base de este estudio. Posteriormente se detallan los módulos de pre-procesado y clasificación que se han combinado en forma de diferentes cadenas de procesamiento orientadas a clasificar datos hiperespectrales de forma supervisada. En este sentido, el núcleo del presente trabajo viene dado por la comparativa de las cadenas de procesamiento consideradas en el marco de dos casos de estudio centrados en la utilización de imágenes hiperespectrales de tejido tumoral, obtenidas mediante sensores hiperespectrales durante intervenciones quirúrgicas y adquiridas de las biopsias cerebrales extraídas durante dichas intervenciones. Como resultado del estudio cuantitativo y comparativo realizado al analizar los resultados de clasificación obtenidos utilizando diferentes cadenas de procesamiento en relación con información de referencia (muestras etiquetadas) disponible para dichas imágenes, se ofrecen una serie de conclusiones y recomendaciones generales acerca del mejor uso posible de los módulos de pre-procesado y clasificación que integran dichas cadenas. Dichas recomendaciones suponen un aspecto innovador en la literatura especializada dedicada al análisis de datos hiperespectrales, y pensamos que serán de gran utilidad para los usuarios de este tipo de datos interesados en aplicaciones relacionadas con la clasificación supervisada de los mismos. Contents 1. Introduction .............................................................................................................................1 1.1. Objectives ..........................................................................................................................1 1.2. Context and motivations ...................................................................................................2 1.3. Petitioner...........................................................................................................................4 1.4. Memory organization ........................................................................................................4 2. Hyperspectral Imaging .............................................................................................................6 2.1. Hyperspectral imaging concept .........................................................................................6 2.2. Medical hyperspectral imaging .........................................................................................9 2.2.1. Brain tumours .............................................................................................................9 2.2.2. Hyperspectral imaging applied to tumour detection ................................................10 2.2.3. Hyperspectral imaging data sets ..............................................................................13 2.2.3.1. Ex vivo samples .................................................................................................14 2.2.3.2. In vivo samples .................................................................................................17 2.3. Hyperspectral analysis techniques ..................................................................................18 2.4. Hyperspectral data classification .....................................................................................20 2.4.1. Unsupervised classification algorithms ....................................................................23 2.4.2. Supervised classification algorithms .........................................................................23 2.4.3. Evaluation metrics used in classification algorithms ................................................24 2.5. Summary .........................................................................................................................29 3. Image processing techniques ................................................................................................30 3.1. Pre-processing techniques for hyperspectral images ......................................................32 3.1.1. De-noising filtering ...................................................................................................32 3.1.1.1. HySime filter .....................................................................................................32 3.1.1.1.1 Noise Estimation ................................................................................34 3.1.1.1.2 Signal Subspace Identification ............................................................35 3.1.1.1.3 HySime Components ..........................................................................37 3.1.1.1.4 Inverse HySime for Hyperspectral image restoration ........................38 3.1.1.2. Smooth filter .....................................................................................................38 3.1.1.2.1 Smoothing algorithms ........................................................................39 3.1.1.2.2 Noise reduction ..................................................................................40 3.1.1.2.3 End effects and the lost points problem ............................................41 3.1.2. Contour detection ....................................................................................................41 3.1.3. Feature extraction techniques..................................................................................48 3.1.3.1. Dimensionality reduction techniques ...............................................................51 3.1.3.1.1. Principal Component Analysis (PCA) .................................................51 3.1.3.1.2. Minimum Noise Fraction (MNF) ........................................................53 3.1.3.1.3. Independent Component Analysis (ICA) ...........................................55 3.1.3.1.4. Locality Preserving Projections (LPP) ................................................60 3.1.3.1.5. Neighborhood Preserving Embedding (NPE) .....................................63 3.1.3.1.6. Kernel Principal Component Analysis (KPCA) ....................................66 3.1.3.2. Spectral unmixing Techniques ..........................................................................71 3.1.3.2.1. Number of endmember estimation ..................................................74 3.1.3.2.2. Endmember extraction algorithms ...................................................75 3.1.3.2.2.1. Orthogonal Subspace Projection (OSP) .............................75 3.1.3.2.2.2. Vertex Component Analysis (VCA) .....................................76 3.1.3.2.2.3. N-FINDR .............................................................................77 3.1.3.2.3. Abundance estimation algorithms ....................................................79 3.1.3.2.3.1. Linear Spectral Unmixing ...................................................79 3.1.3.2.4 Abundance Maps ...............................................................................80 3.1.4. Feature selection ......................................................................................................80 3.1.4.1. Hyperspectral Subspace Identification by Minimum Error ...............................82 3.1.4.2. Virtual Dimensionality (VD) ..............................................................................83 3.2. Classification algorithms ..................................................................................................84 3.2.1. Support Vector Machine Classifier (SVM) .................................................................85 3.2.1.1 Theoretical foundations of SVM classifier .........................................................86 3.2.1.2 Kernel functions available ..................................................................................90 3.3. Summary .........................................................................................................................92 4. Proposed process for cancer detection using hyperspectral images ....................................93 4.1. Process #1: Support vector machine classifier used in brain cancer detection ..............95 4.2. Process #2: Dimensionality reduction techniques used in brain cancer detection ..........97 4.3. Process #3: Spectral unmixing concepts used in brain cancer detection ........................98 4.4. Process #4: Mixed techniques Unmixing-SVM used in brain cancer detection ............100 4.5. Process #5: Different methods for extracting training set using SVM in brain cancer detection ..............................................................................................................................101 4.6. Process #6: Patient simulation using SVM in brain cancer detection ...........................102 4.7. Process #7: Consolidation of processing chains tested using in vivo data in brain cancer detection ..............................................................................................................................103 4.8. Summary .......................................................................................................................105 5. Experimental results ............................................................................................................106 5.1. Results of processing chains ..........................................................................................108 5.1.1. Results of support vector machine classifier used in brain cancer detection .........108 5.1.2. Results of dimensionality reduction techniques used in brain cancer detection ...112 5.1.3. Results of Spectral unmixing concepts used in brain cancer detection ..................131 5.1.4. Results of Mixed techniques Unmixing-SVM used in brain cancer detection .........142 5.1.5. Results of using Different methods for extracting training set using SVM in brain cancer detection...............................................................................................................149 5.1.6. Results of patient simulation using SVM in brain cancer detection........................166 5.1.7. Results of Consolidation of processing chains tested using in vivo data in brain cancer detection...............................................................................................................168 5.2. Global discussion of results ...........................................................................................201 5.3. Summary .......................................................................................................................203 6. Conclusions and future research lines ................................................................................205 6.1. Conclusions ...................................................................................................................205 6.2. Future research lines .....................................................................................................206 Bibliography .............................................................................................................................207 Budget......................................................................................................................................221 Specifications ...........................................................................................................................229 MEMORIA 6 Chapter 2 2. Hyperspectral Imaging 2.1. Hyperspectral imaging concept The observation of a particular object is based on the capture of electromagnetic radiation from the interaction between the object and the radiation source by a measuring instrument or sensor. Electromagnetic radiation received several names depending on the wavelength that characterizes it, as shown in Figure 2.1. For measuring the emitted or reflected radiation by a given surface is necessary to quantify the amount of energy flux that proceeds from the same. For this the radiance measure is used, which depends on factors such as the perceived brightness, reflectance, viewing angles, among others. [5] Spectral detection techniques are based on the fact that all materials in the real world reflect, absorb and emit electromagnetic energy differently in different wavelengths. [6] 7 Figure 2.1: Electromagnetic spectrum. Currently, there is a wide range of instruments or sensors capable of measuring spectral singularities in different wavelengths [7]. The availability of these instruments has provided a redefinition of the digital image concept through the extension of the idea of pixel. It is important to remember that the associated value with each pixel is defined by a numerical value called digital level. This is a numeric value, not visual, but can easily be translated into a visual intensity or gray level by any digital-analog converter. Thus, in a purely spatial schema one pixel is constituted by a single discrete value, while in a spectral schema one pixel consists of a set of values. These values may be understood as Ndimensional vectors [8], where N is the number of spectral bands in which the sensor measures information. Extending the concept of pixel results in what is known as multidimensional image. The order of magnitude of N allows a distinction when talking about multidimensional images. Thus, when the value of N is small, typically a few spectral bands [9], one speaks of multispectral images, whereas when the order of N is hundreds of bands [10], there is talk of hyperspectral imaging, which are also known as hypercubes. 8 Hyperspectral imaging (HSI), like other spectral imaging, collects and processes information from across the electromagnetic spectrum. The goal of hyperspectral imaging is to obtain the spectrum for each pixel in the image of a scene, with the purpose of finding objects, identifying materials, or detecting processes. [11][12] Much as the human eye sees visible light in three bands (red, green, and blue), spectral imaging divides the spectrum into many more bands. This technique of dividing images into bands can be extended beyond the visible as shown in figure [13]. In hyperspectral imaging, the recorded spectra have fine wavelength resolution and cover a wide range of wavelengths. Figure 2.2: Comparison between hypercube and RGB image. Hypercube is a three-dimensional dataset of a two-dimensional image on each wavelength. The lower left is the reflectance curve (spectral signature) of a pixel in the image. RGB color image only has three image bands on red, green, and blue wavelengths respectively. The lower right is the intensity curve of a pixel in the RGB image. Hyperspectral sensors look at objects using a vast portion of the electromagnetic spectrum. Certain objects leave unique fingerprints in the electromagnetic spectrum. If the difference between each of these wavelengths is a few nanometers, it can be obtained for each pixel an “almost continuous” spectral response of the photographed material at that pixel. This response is called spectral signature which is unique to each material. Known as spectral signatures, these fingerprints enable identification of the materials that make up a scanned object and it can be used to identify specific substances in an image. For example, a spectral signature for tumor tissue helps surgeons find tumor areas. Most hyperspectral image pixels are mixed pixels. This is because the spatial resolution of the sensor cannot separate different materials in a pixel. The denotation of mixed pixel is given by the presence of pixels composed by the combination of several pure spectral signatures. In the process of spectral unmixing we find an efficient possibility 9 to express the composition of each mixed pixel in proportion to the pure pixels present. [14][15] With the advantage of acquiring two-dimensional images across a wide range of electromagnetic spectrum, HSI has been applied to numerous areas, engineers have built hyperspectral sensors and processing systems for applications in astronomy, agriculture, biomedical imaging, geosciences, physics, and surveillance. [13] The primary disadvantages of these techniques are cost and complexity. Fast computers, sensitive detectors, and large data storage capacities are needed for analyzing hyperspectral data. Significant data storage capacity is necessary since hyperspectral cubes are large, multidimensional datasets, potentially exceeding hundreds of megabytes. All of these factors greatly increase the cost of acquiring and processing hyperspectral data. 2.2. Medical hyperspectral imaging 2.2.1. Brain tumours Due to the increase in the incidence and mortality from brain tumor in world population in recent decades, the number of research papers related to its diagnosis has grown exponentially. An early diagnosis of such diseases can be vital to prevent a benign tumor evolves into a more aggressive cancer. Brain tumors are due to abnormal growth of cells derived from brain components in the case of primary tumors or tumor cells localized elsewhere in the organism in the case of metastases. According to evolution, size and dimension can be grouped into four degrees of danger. This classification is ranging from grade I tumors (benign tumors of slow growth) to grade IV tumors (malignant cancers with very rapid growth). [16] The 77% of malignant brain cancers belongs to the group of tumors called gliomas. These tumors affect glial cells responsible for supporting neurons and information brain processing. Depending on the type of glial cells affected (astrocytes, oligodendrocytes, ependymal cells, etc.) are different types of gliomas, being multiform glioblastoma (astrocytoma grade IV) the most aggressive and less likely to survive [17][18]. Figure 2.3 shows magnetic resonance imaging of a patient afflicted with a multiform glioblastoma. 10 Figure 2.3: Magnetic resonance of patient suffering from multiform glioblastoma. Symptoms and treatment of these tumors depend largely on the patient's age, the type of tumor and its location in the brain. However, these tumors tend to infiltrate into healthy brain tissue, so that surgery is complex and sometimes impossible. That is why it is vitally important research and study of non-invasive techniques that provide the most accurate diagnosis possible. Proper removal (resection) of brain tumor eliminates malignant tissue and prevents the tumor is reproduced. Increase resection area with an additional safety margin maximizes the chances of patient survival, but in the case of brain tumors, excessive resection can lead to significant damage of all kinds (motor, cognitive, visual, etc.). For this reason, a technique to determine accurately and minimally invasive glioma limits is necessary. 2.2.2. Hyperspectral imaging applied to tumour detection Spectral imaging is a technology that integrates conventional imaging and spectroscopy methods to obtain both spatial and spectral information from an object. Although this technology was originally developed for remote sensing, it has been extended to the biomedical engineering field as a powerful analytical tool for biological and biomedical research [19]. Hyperspectral imaging field is an emerging imaging modality for medical applications. It offers great potential for noninvasive disease diagnosis and surgical 11 guidance. Light delivered to biological tissue undergoes multiple scattering from inhomogeneity of biological structures and absorption primarily in hemoglobin, melanin, and water as it propagates through the tissue [20][21]. One of the important advantages of this technique is that it can acquire reflectance, absorption, or fluorescence spectrum for each pixel in the image. It is assumed that the absorption, fluorescence, and scattering characteristics of tissue change during the progression of disease [22], which can be used to detect the biochemical changes of objects that cannot be identified with traditional gray or color imaging methods [19]. Therefore, the reflected, fluorescent, and transmitted light from tissue captured by HSI carries quantitative diagnostic information about tissue pathology [22-25]. In recent years, advances in hyperspectral cameras, image analysis methods, and computational power make it possible for many exciting applications in the medical field. These applications mainly cover the ultraviolet (UV), visible (VIS), and near-infrared (nearIR or NIR) regions. Table 2.1 defines the spectral range from UV to mid-IR (200 to 25,000 nm) [26]. Short name Full name Spectral range (nm) UV Ultraviolet 200 to 400 VIS Visible 400 to 780 NIR/near-IR Near-infrared 780 to 2500 MIR/mid-IR Mid-infrared 2500 to 25.000 Table 2.1: Spectral range definitions. Visible light penetrates only 1 to 2 mm below the skin and thus obtains information from the sub-papillary [27], while light in the NIR region penetrates deeper into the tissue than VIS or mid-IR radiation [28]. NIR light is preferred for surgical guidance due to its deep penetration into the tissue, which can help the surgeon see through connective tissue for visualizing critical anatomical structures of interest that are not visible and detecting molecules with detectible spectra [13]. By expanding light beyond the visual spectrum, additional information can be obtained to further characterize the cells of interest [29]. HSI acquires a three-dimensional dataset called hypercube, with two spatial dimensions and one spectral dimension. Figure 2.4 [19] shows the concept of the hypercube data captured by a spectral imaging system. These dates can be visualized as a 12 three-dimensional (3-D) cube or a stack of multiple two-dimensional (2-D) images because of its intrinsic structure, in which the cube face is a function of the spatial coordinates and the depth is a function of wavelength [19]. Spatially resolved spectral imaging obtained by HSI provides diagnostic information about the tissue physiology, morphology, and composition [13]. Figure 2.4: The concept of spectral data cube. The data cube contains two spatial dimensions (x and y) and one spectral dimension, in which the cube face is a function of the spatial coordinates and the depth is a function of wavelength. According to the electromagnetic theory, different biochemical constituents commonly have different spectral signatures [30]. These signatures are usually generated by the interactions between materials and electromagnetic waves, such as electron transition, atomic and molecular vibration or rotation. The biological and pathological changes in tissues and organs also have a close relationship with the spectra. Spectral characteristics in different wavelength regions yield a distinguishable spectral signature, making pathological changes distinguishable. Therefore, the spectral imaging technology also can be extended to the biomedical engineering field to estimate the physiological status of biological tissues, since it can take advantage of the spatial relationships among the different spectra in a neighborhood. This technology opens new prospects for life science by which scientists can identify and quantify the relationships among biologically active molecules, observe living organisms noninvasively, perform histopathological and fluorescent analyses, and enhance biological understanding of diseases [19]. 13 Image analysis enables the extraction of diagnostically useful information from a large medical hyperspectral dataset at the tissue, cellular, and molecular levels and is, therefore, critical for disease screening, diagnosis, and treatment. Hypercube with high spatial and spectral resolution may potentially contain more diagnostic information. However, high spatial and spectral dimensions also make it difficult to perform automatic analysis of hyperspectral data. In particular, it is complex in many aspects: (1) high data redundancy due to high correlation in the adjacent bands, (2) variability of hyperspectral signatures, and (3) curse of dimensionality [7]. With abundant spatial and spectral information available, advanced image classification methods for hyperspectral datasets are required to extract, unmix, and classify relevant spectral information. The goal is not only to discriminate between different tissues (such as healthy and malignant tissue) and provide diagnostic maps, but also to decompose mixtures into the spectra of pure molecular constituents and correlate these molecular fingerprints (biomarkers) with disease states. Although hyperspectral image analysis methods have been intensively investigated in the remote sensing area, their development and application in medical domain lag far behind. The relationships between spectral features and underlying biomedical mechanisms are not well understood. The basic steps for hyperspectral image analysis generally involve preprocessing, feature extraction and feature selection, and unmixing and/or classification [13]. 2.2.3. Hyperspectral imaging data sets Throughout this dissertation two different data are used: ex vivo samples and in vivo samples. In the first stage of this work, all the experiments were performed with ex vivo samples. In the second stage, once the in vivo samples were obtained, the main processing chains were repeated using such data. This section explains how datasets used along this work were obtained. 14 2.2.3.1. Ex vivo samples The ex vivo samples are formed by living tissues which were removed from a real patient in a surgery. Then, several pictures were taken with a hyperspectral camera from different angles. The camera used has 1040 spectral bands ranging from 380 to 1028 nanometers. For each picture the hyperspectral data cube is formed from the files generated by the hyperspectral camera. In this regard, we have 11 hyperspectral data cube, which have 3 different angles and can contain up to 3 different tissue types: healthy, tumor and necrosis. The process to obtain the samples is detailed below: 1. The hyperspectral data cube is formed from the files generated by the hyperspectral camera. Figure 2.5: Hypercube 2 band 350. 2. The image is cropped in order to it can be processed more easily thanks to its reduced size and for separating the different samples (healthy, tumor and necrosis) in each image. 15 Figure 2.6: Hypercube 2 band 350 separated by classes: healthy, tumor and necrosis. 3. The samples are filtered because these contain a lot of noise. The filters used are Hyperspectral Subspace Identification by Minimum Error (HySime) and smooth, and these are explained in detail in the following sections. 4. Finally the contour of the sample is determined by removing the rest of the image, in order to determine the number of samples of each type at our disposal. The contour functions are explained in detail in the following sections. Figure 2.7: Tumor sample hypercube 2 without filter, HySime filter and smooth filter. 22  In the targets detection the pixels are classified into two classes called “target” and “background”, depending on whether they contain the target sought or not. The binary classification problem can be formulated mathematically as follows. Let R the N-dimensional space formed by all pixels in the hyperspectral image. Given a Ndimensional array 𝒖=(𝑢1,𝑢2,…,𝑢𝑁)𝑇 associated with a specific pixel. The binary classification consists of dividing the space R into two regions, 𝑅𝑜 and 𝑅𝑓, such that 𝒖 is classified as “target” if 𝒖∈𝑅𝑜 and “background” if 𝒖∈𝑅𝑓. This problem can be illustrated graphically using a scatterplot between two bit correlated bands of hyperspectral imaging, as shown in Figure 2.10. Figure 2.10: Graphic illustration of the classification problem in hyperspectral imaging. As shown in Figure 2.10, the ideal situation in a classification problem occurs when the separation between target and background is clearly defined in clearly distinguishable clusters. 23 The hyperspectral image classification algorithms can be divided into two broad categories [37]: supervised algorithms and unsupervised algorithms:  Unsupervised algorithms: These algorithms assume that there is no a priori knowledge about existing classes. The aim of these techniques is to identify automatically classes or groups of pixels using a similarity metric.  Supervised algorithms: It start with some knowledge about existing classes, from which it can be derived some classification criteria. This approach tends to be given by a preliminary step in which spectral signatures of the existing classes are selected. 2.4.1. Unsupervised classification algorithms The pixels classification techniques in an unsupervised way for hyperspectral imaging are in full development phase [7]. Among the existing techniques highlights the KMeans method [40], which supposes the existence of K classes (parameter to be determined a priori) and performs a grouping of the pixels of the image in such classes, using purely statistical methods of vector quantization based on the average spectra of these classes. Moreover, the ISODATA method [41] also requires initialization of a K parameter relating to the number of desired classes, prior to the execution of the algorithm. In addition, this method requires information on the minimum number of signatures belonging to one class. If the initial value of K is low, the dispersion between classes can be very high. Conversely, if the initial value of K is high, the distance between classes can be very small, causing the partitioning of the same class in several classes similar to each other. Overall, the recent literature shows that the results obtained by these two techniques have not been very satisfactory, except in very specific applications [42]. 2.4.2. Supervised classification algorithms In supervised classification techniques highlights the matched filters and Spectral Angle Mapper (SAM) method, both based on first-order statistics [43]. Within this category may also be other classifiers as the nearest neighbor, minimum distance, parallelepiped or maximum likelihood techniques. However, this dissertation is focus on the Support Vector 24 Machine (SVM) technique which has demonstrated excellent performance when working with high-dimensional data such as hyperspectral data. Once presented the most common techniques of hyperspectral image classification, this section is concluded highlighting some techniques used to evaluate the performance of these algorithms. 2.4.3. Evaluation metrics used in classification algorithms The large number of existing techniques as well as the continued proliferation of new methodologies makes clear the need for comparative schemes or metrics to qualitatively analyze the performance of new methodologies rose, contrasting its results with those provided by the existing ones. Most of the evaluation techniques of digital image analysis algorithms are based on the concept of ground-truth [44]. Ideally it is possible to define the concept of ground-truth as the optimal classification or interpretation result which should get an algorithm [45][46]. The ground-truth usually comes characterized by relevant information about the properties in the real world of a set of targets that are desired to identify or characterize. This information is usually obtained by measurements made directly in the study area covered by the image [47], although it is also possible to obtain ground-truth information by applying algorithmic techniques [6]. In any case, the first alternative is the more reliable, but it can be expensive due to the need to label all the samples obtained [48]. Assuming the existence of ground-truth information, there are several methodologies to compare that information with the results provided by an image analysis algorithm. This section provides a brief description of the different metrics that can be applied to assess the ability of a computer algorithm in terms of classification and identification of targets of interest in a digital image. In particular, it will highlight one of the most widely used approaches called the confusion matrix, that will serve for the subsequent comparative study of this dissertation and from which other metrics are derived as the percentage of success in the classification. 25 The confusion matrix [49] is a technique to assess the accuracy of classification algorithms of digital images. This technique assumes that the ground-truth information is expressed as a thematic map [50][51] characterized by the following properties: 1. Each pixel is labeled as belonging to a particular class, so as to have N classes or reference regions {𝑅𝑖}𝑖=1 𝑁. 2. The reference regions are mutually exclusive of each other, ie, two different regions have no pixel in common: 𝑅𝑖∩𝑅𝑗=∅,∀𝑖≠𝑗. Assuming that each pixel i of the image to evaluate I is assigned by the algorithm as belonging to a certain class 𝐶𝑖, so as to have N classes. 𝐶𝑖 sets are a partition of the image to be evaluated, that is, the union of all sets resulting in the image and two different sets have no element in common: ⋃𝐶𝑖=𝐼 𝑎𝑛𝑑 𝐶𝑖∩𝐶𝑗=∅,∀𝑖≠𝑗 𝑁 𝑖=1 (2.1) Given the above considerations, Figure 2.11 shows an example of the process of constructing a confusion matrix. In the figure, the thematic map of ground-truth classification associated with the image to classify, the classification result provided by a given algorithm to that image, and the confusion matrix that quantifies the accuracy of the algorithm in the classification task are shown. As shown, the matrix entries are expressed in the form of 𝑎𝑗𝑘 where 𝑎𝑗𝑘=𝑐𝑎𝑟𝑑𝑖𝑛𝑎𝑙{𝐶𝑗∩𝑅𝑘}, that is, the number of pixels of the resulting region when performing the intersection between a 𝐶𝑗 class obtained by the algorithm and a ground-truth class 𝑅𝑘 [49]. 26 Figure 2.11: Example of constructing of a confusion matrix. From the confusion matrix can be derived some generic accuracy measures [52] as:  Overall accuracy (OA): Percentage of pixels correctly classified in all classes: 𝑂𝐴=∑𝑎𝑖𝑖 𝑁 𝑖 ∑𝑎𝑖𝑗 𝑁 𝑖𝑗 ×100 (2.2)  Average accuracy (AA): Means percentage of classification accuracy per class for all classes: 𝐴𝐴=𝑎𝑖𝑖 ∑𝑎𝐴𝑖 𝑁 𝑗×100 (2.3) AA can also be calculated only for a particular class instead of as the mean of all the classes. In that case it is known as the Average accuracy per class. OA and AA are percentages should approach 100%, which would be the ideal classification. When the reference set is not well defined, the OA will not be representative regarding the true performance of the classifier. For example, if a class has very few pixels of reference, its influence on the computation of OA will be very low, while the AA will become more important, since it is an average done with the number of classes not with the total number of pixels. If the differences arising between these two measures 27 are high, then it may indicate that there has been a poor classification for a specific class, which can affect the overall results of classification. Moreover, there are two additional metrics really important in the medical field: sensitivity and specificity. Let us imagine a study evaluating a new test that screens healthy and tumor samples. Each sample taking the test either has or does not have the disease. The test outcome can be positive (predicting that the sample is a tumor sample) or negative (predicting that the sample is not a tumor sample and therefore it is healthy sample). The test results for each sample may or may not match the sample's actual status. In that setting:  True positive: Tumor sample correctly diagnosed as tumor.  False positive: Healthy sample incorrectly identified as tumor.  True negative: Healthy sample correctly identified as healthy.  False negative: Tumor sample incorrectly identified as healthy. In general, Positive = identified and negative = rejected. Therefore:  True positive = correctly identified.  False positive = incorrectly identified.  True negative = correctly rejected.  False negative = incorrectly rejected. The four outcomes can be formulated in a confusion matrix as follows: Condition (ground-truth) Total population Condition positive Condition negative Test outcome positive True positive False positive (Type I error) Test outcome negative False negative (Type II error) True negative Table 2.4: Confusion matrix. Sensitivity and specificity are statistical measures of the performance of a binary classification test: 28  Sensitivity (true positive rate): Relates to the test's ability to correctly detect patients who do have a condition. Consider the example of a medical test used to identify a disease. Sensitivity of the test is the proportion of samples known to have the disease, which test positive for it. Mathematically, this can be expressed as: 𝑆𝑒𝑛𝑠𝑖𝑡𝑖𝑣𝑦= 𝑇𝑟𝑢𝑒 𝑃𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑠 𝑇𝑟𝑢𝑒 𝑃𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑠 + 𝐹𝑎𝑙𝑠𝑒 𝑁𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑠= =𝑇𝑟𝑢𝑒 𝑃𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑠 𝑇𝑜𝑡𝑎𝑙 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 𝑖𝑛 𝑡ℎ𝑒 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑜𝑟 𝑖𝑛 𝑡ℎ𝑒 𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ =𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 𝑜𝑓 𝑎 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 𝑡𝑒𝑠𝑡 𝑔𝑖𝑣𝑒𝑛 𝑡ℎ𝑎𝑡 𝑡ℎ𝑒 𝑝𝑎𝑡𝑖𝑒𝑛𝑡 𝑖𝑠 𝑖𝑙𝑙 (2.4) A test with high sensitivity has a low type II error rate.  Specificity (true negative rate): Relates to the test's ability to correctly detect patients without a condition. Consider the example of a medical test for diagnosing a disease. Specificity of a test is the proportion of healthy samples known not to have the disease, which will test negative for it. Mathematically, this can also be written as: 𝑆𝑝𝑒𝑐𝑖𝑓𝑖𝑐𝑖𝑡𝑦= 𝑇𝑟𝑢𝑒 𝑁𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑠 𝑇𝑟𝑢𝑒 𝑁𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑠 + 𝐹𝑎𝑙𝑠𝑒 𝑃𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑠= =𝑇𝑟𝑢𝑒 𝑁𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑠 𝑇𝑜𝑡𝑎𝑙 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 𝑖𝑛 𝑡ℎ𝑒 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑜𝑟 𝑖𝑛 𝑡ℎ𝑒 𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ =𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 𝑜𝑓 𝑎 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 𝑡𝑒𝑠𝑡 𝑔𝑖𝑣𝑒𝑛 𝑡ℎ𝑎𝑡 𝑡ℎ𝑒 𝑝𝑎𝑡𝑖𝑒𝑛𝑡 𝑖𝑠 𝑤𝑒𝑙𝑙 (2.5) A test with a high specificity has a low type I error rate. For any test, there is usually a trade-off between the measures. A perfect predictor would be described as 100% sensitive (e.g., all sick are identified as sick) and 100% specific (e.g., all healthy are not identified as sick); however, theoretically any predictor will possess a minimum error bound known as the Bayes error rate. Finally, it is defined the Kappa coefficient (K), that is a statistical measure to assess the reliability of an agreement between a fixed number of evaluators to assign categorical classifications to a number of elements which are classified by these evaluators. This measure calculates the degree of agreement in the classification what would be expected by random effect, i.e., when evaluators are not absolutely sure and just venture an answer. 29 𝑘=𝑐𝑜𝑟𝑟𝑒𝑐𝑡 𝑐𝑙𝑎𝑠𝑠𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 − 𝑎𝑔𝑟𝑒𝑒𝑚𝑒𝑛𝑡 𝑏𝑦 𝑐ℎ𝑎𝑛𝑐𝑒 𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 1 − 𝑎𝑔𝑟𝑒𝑒𝑚𝑒𝑛𝑡 𝑏𝑦 𝑐ℎ𝑎𝑛𝑐𝑒 𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 (2.6) The kappa coefficient can be interpreted considering the following table: Kappa coefficient (K) Classification can be regarded as Below 0.4 Poor 0.41-0.60 Moderate 0.61-0.75 Good 0.76-0.80 Excellent 0.81 and above Almost perfect Table 2.5: Interpretation of Kappa coefficient. If evaluators are in full agreement, then 𝑘=1. If no agreement among evaluators (other than what would be expected by chance) then 𝑘=0 [53]. There are other metrics such as the commission errors and omission errors that are not considered in this study. 2.5. Summary Hyperspectral imaging is a very powerful tool that can be used to provide a solution in many different fields. In this dissertation it is intended to apply hyperspectral imaging in the medical field, particularly in the brain cancer detection. For this purpose, a study of main existing techniques for hyperspectral imaging preprocessing is performed because it is necessary to properly process hyperspectral data with the ultimate goal of obtaining the best possible classification. In order to evaluate the classification process, a set of standard evaluation metrics will be used to subsequently to draw conclusions with regard to the efficiency approaches for brain cancer detection. 30 Chapter 3 Image processing techniques Image processing is the processing of images using mathematical operations in any form of signal processing for which the input is an image, such as a photograph, video frame or hyperspectral imaging; the output of image processing may be either an image or a set of characteristics or parameters related to the image [54]. Namely, image processing is the set of techniques applied to images in order to improve quality or facilitate the search for information. The processing chain in this project can be divided into five different stages, as we can observe in Figure 3.1: Figure 3.1: Stages of hyperspectral data processing chain followed in this project. The classification system is composed by the preprocessing stage and the classification algorithms. Firstly, the initial set of samples is studied in order to delete the anomalous samples taken during the surgery. In this stage, the hyperspectral cubes are spatially reduced in size, in order to process the region of interest. 31 On the one hand, in the case of ex vivo samples, its contour is obtained with the purpose of extracting the samples from the rest of the image (table, test plates, shadows and glow in the samples), and to prepare the data for the classification process in response to the inclination of the samples (no inclination, 23.58 degrees and 19.10 degrees) and according to the kind of filter used. For this reason, in this stage we use two types of filter: HySime and Smooth, as well as a function to obtain the contour. On the other hand, in the case of in vivo samples, the images are formed by several elements, labeled directly on the image thereby forming the in vivo dataset. Then, these samples are processed as the ex vivo samples. Features are functions of the original measurement variables that are useful for classification and/or pattern recognition [55]. The following step of the processing chain, feature extraction, is the process of defining a set of features, or image characteristics, which will most efficiently or meaningfully represent the information that is important for analysis and classification. In this stage there are several applicable methodologies, of which will be applied in this dissertation, mainly dimensionality reduction and unmixing techniques. The next stage of the chain relies on the feature selection, also important in reducing the high dimensionality of the data. The fact of choosing a set of features whose dimensionality is the most appropriate and reasonable is a hot topic in the scientific literature, so efficient and fast algorithms that perform the process of combining bands or characteristics for a given problem are needed. This is a complex process that cannot be defined with a unique approach. There are several techniques of statistical indices for feature selection that will be seen in detail in the corresponding section. Finally, once the data have been reduced and/or processed, and with some analyzes already performed, they will be classified. As noted in the previous chapter, classification is generally a process in which the individual elements or items are differentiated into groups, based on quantitative information of one or more inherent characteristics of the elements, usually by items previously labeled using sets or training patterns. Finally, and as an option, the resulting 38 3.1.1.1.4 Inverse HySime for Hyperspectral image restoration Once the HySime components have been obtained and the noise segregated the Hyperion data can be restored to its original spectral space without noise. The noise free original spectral space consisting of the noise less signals only can be achieved by performing an inverse HySime transform. As the original image data Y was transformed into HySime components, 𝑌𝐻𝑦𝑆𝑖𝑚𝑒 in the HySime space, the inverse transformation can be achieved by inverting the projection matrix 𝐸 and multiplying it with𝑌𝐻𝑦𝑆𝑖𝑚𝑒. The expression for restoration of signals to original spectral space is as follows, 𝑍=(𝐸𝑇)−1𝑌𝐻𝑦𝑆𝑖𝑚𝑒 (3.19) where, 𝑌𝐻𝑦𝑆𝑖𝑚𝑒=(𝑦1𝐻𝑦𝑆𝑖𝑚𝑒,……,𝑦𝑝𝐻𝑦𝑆𝑖𝑚𝑒,0,……,0) - is the HySime component image p - the number of 1st components selected for restoring to the original spectral space and 𝑍 - the restored image 3.1.1.2. Smooth filter In many experiments in physical science, the true signal amplitudes (y-axis values) change rather smoothly as a function of the x-axis values, whereas many kinds of noise are seen as rapid, random changes in amplitude from point to point within the signal. In the latter situation it may be useful in some cases to attempt to reduce the noise by a process called smoothing. In smoothing, the data points of a signal are modified so that individual points that are higher than the immediately adjacent points (presumably because of noise) are reduced, and points that are lower than the adjacent points are increased, thus is achieved that the signal is more homogeneous. This naturally leads to a smoother signal. As long as the true underlying signal is actually smooth, then the true signal will not be much distorted by smoothing, but the noise will be reduced [62]. When processing a signal by the Smooth filter, the goal is to create an approximate function that attempts to capture important signal patterns, leaving out the noise. The Smooth filter smooths the signal by using a moving average filter. In the moving average filters each output value is obtained as the average of a subset of the original data. Its purpose is to highlight the significant pattern i.e. it reduce the noise smoothing the fluctuations in short periods, thus highlighting trends or long periods cycles [62]. 39 Moving average filters [63] are the most common filters used in signal processing, mostly because it is the easiest to understand and use digital filter. Despite its simplicity, it is optimal for the task of reducing random noise. This makes them the best filters to signals in the time domain. Instead, moving average filters are the worst for signals encoded in the frequency domain, as they have little ability to separate frequency bands from each other [64]. 3.1.1.2.1 Smoothing algorithms Most smoothing algorithms are based on the "shift and multiply" technique, in which a group of adjacent points in the original data are multiplied point-by-point by a set of numbers (coefficients) that defines the smooth shape, the products are added up to become one point of smoothed data, then the set of coefficients is shifted one point down the original data and the process is repeated. The simplest smoothing algorithm is the rectangular or unweighted sliding-average smooth; it simply replaces each point in the signal with the average of m adjacent points, where m is a positive integer called the smooth width. For example, for a 3-point smooth (m = 3): 𝑆𝑗=𝑌𝑗−1+𝑌𝑗+𝑌𝑗+1 3 (3.20) for j = 2 to n-1, where Sj the jth point in the smoothed signal, Yj the jth point in the original signal, and n is the total number of points in the signal. Similar smooth operations can be constructed for any desired smooth width, m. Usually m is an odd number. If the noise in the data is "white noise" (that is, evenly distributed over all frequencies) and its standard deviation is s, then the standard deviation of the noise remaining in the signal after the first pass of an unweighted sliding-average smooth will be approximately s over the square root of m (s/sqrt(m)), where m is the smooth width. The triangular smooth is like the rectangular smooth, above, except that it implements a weighted smoothing function. For a 5-point smooth (m = 5): 𝑆𝑗=𝑌𝑗−2+2𝑌𝑗−1+3𝑌𝑗+2𝑌𝑗+1+𝑌𝑗+2 9 (3.21) for j = 3 to n-2, and similarly for other smooth widths. 40 It is often useful to apply a smoothing operation more than once, that is, to smooth an already smoothed signal, in order to build longer and more complicated smooths. For example, the 5-point triangular smooth above is equivalent to two passes of a 3-point rectangular smooth. Three passes of a 3-point rectangular smooth result in a 7point "pseudo-Gaussian" or haystack smooth, for which the coefficients are in the ratio 1 3 6 7 6 3 1. The general rule is that n passes of a w-width smooth results in a combined smooth width of n*w-n+1. For example, 3 passes of a 17-point smooth results in a 49point smooth. These multipass smooths are more effective at reducing high-frequency noise in the signal than a rectangular smooth. In all these smooths, the width of the smooth m is chosen to be an odd integer, so that the smooth coefficients are symmetrically balanced around the central point, which is important because it preserves the x-axis position of peaks and other features in the signal. (This is especially critical for analytical and spectroscopic applications because the peak positions are often important measurement objectives). The Savitzky-Golay smooth is based on the least-squares fitting of polynomials to segments of the data. Compared to the sliding-average smooths, the Savitzky-Golay smooth is less effective at reducing noise, but more effective at retaining the shape of the original signal. It is capable of differentiation as well as smoothing. The algorithm is more complex and the computational times are greater than the smooth types discussed above, but with modern computers the difference is not significant [62][65]. 3.1.1.2.2 Noise reduction Smoothing usually reduces the noise in a signal. If the noise is “white” (that is, evenly distributed over all frequencies) and its standard deviation is s, then the standard deviation of the noise remaining in the signal after one pass of a triangular smooth will be approximately s*0.8/sqrt(m), where m is the smooth width. Smoothing operations can be applied more than once: that is, a previously-smoothed signal can be smoothed again. In some cases this can be useful if there is a great deal of high-frequency noise in the signal. However, the noise reduction for white noise is less in each successive smooth. For example, three passes of a rectangular smooth reduces white noise by a factor of approximately s*0.7/sqrt(m), only a slight improvement over two passes. 41 The frequency distribution of noise, designated by noise color, substantially affects the ability of smoothing to reduce noise. Because smoothing is a low-pass filter process, it affects low frequency (pink) noise less, and high-frequency (blue) noise more, than white noise. It should be clear that smoothing can never completely eliminate noise, because most noise is spread out over a wide range of frequencies, and smoothing simply reduces the noise in part of its frequency range. Only for some very specific types of noise (e.g. discrete frequency noise or single-point spikes) is there hope of anything close to complete noise elimination. 3.1.1.2.3 End effects and the lost points problem Note in the equations above that the 3-point rectangular smooth is defined only for j = 2 to n-1. There is not enough data in the signal to define a complete 3-point smooth for the first point in the signal (j = 1) or for the last point (j = n), because there are no data points before the first point or after the last point (similarly, a 5-point smooth is defined only for j = 3 to n-2, and therefore a smooth cannot be calculated for the first two points or for the last two points). In general, for an m-width smooth, there will be (m-1)/2 points at the beginning of the signal and (m-1)/2 points at the end of the signal for which a complete m-width smooth cannot be calculated. There are two approaches for this problem. One is to accept the loss of points and trim off those points or replace them with zeros in the smooth signal. The other approach is to use progressively smaller smooths at the ends of the signal, for example to use 2, 3, 5, 7... points smooths for signal points 1, 2, 3, and 4..., and for points n, n-1, n-2, n-3..., respectively. The later approach may be preferable if the edges of the signal contain critical information, but it increases execution time [65]. 3.1.2. Contour detection Figure 3.4: Stages of hyperspectral data processing chain. This process is based on contour detection of tumor, healthy and necrosis samples. In fact, the contour of the samples is delimited by a function in order to highlight the 42 samples ignoring the rest of the elements in the image (table, test plates and glow in the samples). This is possible due to the difference in density between different tissues and other elements in the image. This function take a representative band which the image can be seen clearly, then it takes each pixel in the image and assigns a value of 1 or 0 depending on whether there is sample or not. To differentiate the samples from the rest of the image a contour function belonging to a MATLAB toolbox is used. This way, we obtain a binary image which has 1 in the samples and 0 in the rest of the image. This function has a condition with a decision threshold to decide if the pixel contains a sample or not. This threshold was chosen by trial and error as shown in the following images. Figure 3.5: Tumor sample 3 without inclination according to the threshold. In figure 3.5, the threshold is too low and we can see how some samples are lost. According as the threshold increases in order to compensate for sample loss, we can see that there are areas without samples which are marked as such. In the end, after several tests, the threshold was established in 40 because this value generally provides good average results. 43 (a) (b) Figure 3.6: Tumor 3 without inclination: (a) band 350 and (b) contour with threshold=40. However, this function is only valid for samples without inclination. Because in the samples with inclination, the samples generate a shadow which is confused with the sample and it produces a false samples extraction how is shown in figure 3.7. (a) (b) Figure 3.7: Tumor 7 19.10 degrees of inclination: (a) band 350 and (b) contour with threshold=40. To solve this problem, it is necessary to show the image in other format which allow to differentiate clearly the sample from the shadow. For this purpose, firstly the image is changed to RGB format and then it is changed to HSV format. Due to these changes in the input data, it is necessary to modify the condition and the decision threshold in the contour function. This condition and its threshold were chosen by trial and error as shown in the images of Figure 3.8. 44 Figure 3.8: Tumor sample 7, 19.10 degrees of inclination according to the threshold. In the first set of images, the threshold chosen is composed by a low value as we can see in the figure 3.8. Several studies were performed in order to find the correct value in which the experimental samples (tumor, necrosis and healthy) could be extracted. In this set of experiments, the threshold was established to 0.67. (a) (b) Figure 3.9: Tumor 7 19.10 degrees of inclination: (a) band 350 HSV and (b) contour with threshold=0.67. Our study, conducted using three kind of inclination collected by different filters, provides several set of samples. In this regard, the nomenclatures used in this dissertation 45 are: inclination 1, inclination 2 and inclination 3 to designate the inclinations of 0, 23.58, 19.10 degrees respectively, class 1, class 2 and class 3 to refer healthy, tumor and necrosis samples respectively, and filter 1 and filter 2 to denote HySime and smooth filter respectively. The results of preprocessing stage are shown in the table 3.2, where we can observe the number of samples for each case: Inclination Filter Samples Class 1 Class 2 Class 3 Total Inclination 1 No-Filter 2917 5486 1530 9933 Filter 1 2971 5543 1555 10069 Filter 2 2986 5258 1403 9647 Inclination 2 No-Filter 811 1787 441 3039 Filter 1 807 1785 443 3035 Filter 2 741 1789 423 2953 Inclination 3 No-Filter 802 3839 1464 6105 Filter 1 803 3827 1469 6099 Filter 2 742 3834 1404 5980 Table 3.1: Samples obtained. Table 3.2 shows the number of samples corresponding to 3 %, 5% and 10% of the global set of samples. They are obtained by randomly selecting of the global pixels. 46 Inclination Filter Training 3% 5% 10% Inclination 1 No-Filter 300 498 996 Filter 1 303 504 1008 Filter 2 291 483 966 Inclination 2 No-Filter 93 153 306 Filter 1 93 153 306 Filter 2 90 150 297 Inclination 3 No-Filter 144 240 480 Filter 1 183 306 612 Filter 2 180 300 600 Table 3.2: Trainings generated. In this set we have 3 possible classes: the number of samples per class for each training set is equal to the total number of samples in such training divided by 3. Moreover, each test is formed by removing the training samples from the total samples, so each test is composed by the number of samples resulting from subtracting the total number of samples for each case minus the number of samples of a particular training. This process is shown in figure 3.10. 47 Figure 3.10: Representation of the generation of samples set, training samples and test samples. In order to study all of the cases, the data with different angles are mixed together. In this case we remove the dependence over the angles. Filter Samples No-Filter 19077 Filter 1 19203 Filter 2 18580 Table 3.3: Total samples. Currently, the resulting trainings for these datasets are obtained. Test samples Samples set Class 1 Training samples Class 2 Class 3 Randomly selected 54 3. Next, a set of components containing weighted information about the variance in the original data set is obtained. For this is used a MNF index that estimates the ratio between signal and noise present in the components provided by the PCA transform. The component having the minimum noise fraction is that whose associated eigenvector, 𝒗, maximizes the expression (3.24). The main difference between the PCA transformation and MNF transformation is the fact that in the second case a more detailed description of the relationship between the amount of signal and the amount of noise present in the image is performed [80]. Thus, the first band resulting from the MNF transformation is the one with highest SNR ratio. The second band has a better SNR than the third, and so on. Figure 3.14: Application example of the MNF transform on a hyperspectral image. As a consequence of the more accurate estimation of the noise conditions present in the image, in certain applications the MNF decomposition may provide more robust results than the PCA transform [81], because it is less sensitive to outliers and noisy pixels. The MNF transform performs a translation of the data, so that the coordinate origin is the 55 centroid of the resulting point cloud. Sometimes, this feature allows obtaining a better description of the data. By comparison, Figure 3.14 [5] shows the first 20 bands obtained from the application of the MNF transform on the real hyperspectral image previously used. 3.1.3.1.3. Independent Component Analysis (ICA) In signal processing, Independent Component Analysis (ICA) is a statistical and computational method for separating a multivariate signal into additive subcomponents, revealing hidden factors that underlie sets of random variables, measurements, or signals. This is done by assuming that the subcomponents are non-Gaussian signals and that they are statistically independent from each other. Independent Component Analysis attempts to decompose a multivariate signal into independent non-Gaussian signals. As an example, sound is usually a signal that is composed of the numerical addition, at each time t, of signals from several sources. The question then is whether it is possible to separate these contributing sources from the observed total signal. When the statistical independence assumption is correct, blind ICA separation of a mixed signal gives very good results. It is also used for signals that are not supposed to be generated by a mixing for analysis purposes. An important note to consider is that if N sources are present, at least N observations are needed to recover the original signals. This constitutes the square case (J = D, where D is the input dimension of the data and J is the dimension of the model). Other cases of underdetermined (J > D) and overdetermined (J < D) have been investigated. The fact that the ICA separation of mixed signals provides very good results are based on two assumptions and three effects of mixing source signals. The two assumptions are: 1. The source signals are independent of each other. 2. The values in each source signal have non-Gaussian distributions. The three effects of mixing source signals are: 1. Independence: As per assumption 1, the source signals are independent; however, their signal mixtures are not. This is because the signal mixtures share the same source signals. 56 2. Normality: According to the Central Limit Theorem, the distribution of a sum of independent random variables with finite variance tends towards a Gaussian distribution. Loosely speaking, a sum of two independent random variables usually has a distribution that is closer to Gaussian than any of the two original variables. Here we consider the value of each signal as the random variable. 3. Complexity: The temporal complexity of any signal mixture is greater than that of its simplest constituent source signal. Those principles contribute to the basic establishment of ICA. If the signals we happen to extract from a set of mixtures are independent like sources signals, or have non-Gaussian histograms like source signals, or have low complexity like source signals, then they must be source signals. [82][83] ICA defines a generative model for the observed multivariate data, which is typically given as a large database of samples. In the model, the data variables are assumed to be linear mixtures of some unknown latent variables, and the mixing system is also unknown. The latent variables are assumed to be non-Gaussian and mutually independent and they are called the independent components of the observed data. These independent components, also called sources or factors, can be found by ICA. ICA finds the independent components (also called factors, latent variables or sources) by maximizing the statistical independence of the estimated components. We may choose one of many ways to define independence, and this choice governs the form of the ICA algorithm. The two broadest definitions of independence for ICA are: 1. Minimization of mutual information (uses measures like Kullback-Leibler Divergence and maximum entropy). 2. Maximization of non-Gaussianity (motivated by the central limit theorem, uses kurtosis and negentropy). Typical algorithms for ICA use centering (subtract the mean to create a zero mean signal), whitening (usually with the eigenvalue decomposition), and dimensionality reduction as preprocessing steps in order to simplify and reduce the complexity of the problem for the actual iterative algorithm. Whitening and dimension reduction can be achieved with Principal Component Analysis (PCA). Whitening ensures that all dimensions are treated equally a priori before the algorithm is run. Well-known algorithms for ICA 57 include infomax, FastICA, and JADE, but there are many others. In this case, ICA JADE was used for this project. Linear independent component analysis can be divided into noiseless and noisy cases, where noiseless ICA is a special case of noisy ICA. The data are represented by the random vector: 𝑥=(𝑥1,…,𝑥𝑚)𝑇 (3.25) and the components as the random vector: 𝑠=(𝑠1,…,𝑠𝑛)𝑇 (3.26) The task is to transform the observed data 𝑥, using a linear static transformation 𝑊 as: 𝑠=𝑊𝑥 (3.27) into maximally independent components 𝑠 measured by some function: 𝐹(𝑠1,…,𝑠𝑛) (3.28) of independence. In the linear noiseless ICA, the components 𝑥𝑖 of the observed random vector: 𝑥=(𝑥1,…,𝑥𝑚)𝑇 (3.29) are generated as a sum of the independent components 𝑠𝑘,𝑘=1,…,𝑛: 𝑥𝑖=𝑎𝑖,1𝑠1+⋯+𝑎𝑖,𝑘𝑠𝑘+⋯+𝑎𝑖,𝑛𝑠𝑛 (3.30) weighted by the mixing weights 𝑎𝑖,𝑘. The same generative model can be written in vectorial form as: 𝑥=∑𝑠𝑘𝑎𝑘 𝑛 𝑘=1 (3.31) where the observed random vector 𝑥 is represented by the basis vectors: 𝑎𝑘=(𝑎1,𝑘,…,𝑎𝑚,𝑘)𝑇 (3.32) 58 The basis vectors 𝑎𝑘 form the columns of the mixing matrix: 𝐴=(𝑎1,…,𝑎𝑛) (3.33) and the generative formula can be written as: 𝑥=𝐴𝑠 (3.34) where 𝑠=(𝑠1,…,𝑠𝑛)𝑇 (3.35) Given the model and realizations (samples) 𝑥1,…,𝑥𝑁 of the random vector 𝑥, the task is to estimate both the mixing matrix 𝐴 and the sources 𝑠. This is done by adaptively calculating the 𝑤 vectors and setting up a cost function which either maximizes the nongaussianity of the calculated: 𝑠𝑘=(𝑤𝑇∗ 𝑥) (3.36) or minimizes the mutual information. In some cases, a priori knowledge of the probability distributions of the sources can be used in the cost function. The original sources 𝑠 can be recovered by multiplying the observed signals 𝑥 with the inverse of the mixing matrix: 𝑊=𝐴−1 (3.37) also known as the unmixing matrix. Here it is assumed that the mixing matrix is square (𝑛=𝑚). If the number of basis vectors is greater than the dimensionality of the observed vectors, 𝑛>𝑚, the task is overcomplete but is still solvable with the pseudo inverse. In the case of linear noisy ICA with the added assumption of zero-mean and uncorrelated Gaussian noise: 𝑛~𝑁(0,𝑑𝑖𝑎𝑔(𝛴)) (3.38) the ICA model takes the form: 𝑥=𝐴𝑠+ 𝑛 (3.39) The ICA transformation is illustrated graphically in Figure 3.15 [84]. 59 Figure 3.15: Graphic illustration of the ICA transformation. An example of the application of the ICA transformation to real hyperspectral image is shown in Figure 3.16. 60 Figure 3.16: Application example of the ICA transform on a hyperspectral image. ICA is somehow related to Principal Component Analysis (PCA). What distinguishes ICA from other methods is that it looks for components that are both statistically independent, and non-Gaussian. ICA is a much more powerful technique, however, capable of finding the underlying factors or sources when these classic methods fail completely [85]. 3.1.3.1.4. Locality Preserving Projections (LPP) Locality Preserving Projections (LPP) are linear projective maps that arise by solving a variational problem that optimally preserves the neighborhood structure of the data set. LPP should be seen as an alternative to classical linear technique based on Principal Component Analysis (PCA). When the high dimensional data lies on a low dimensional manifold embedded in the ambient space, the Locality Preserving Projections are obtained by finding the optimal linear approximations to the eigenfunctions of the Laplace Beltrami operator on the manifold. As a result, LPP shares many of the data representation 61 properties of nonlinear techniques such as Laplacian Eigenmaps or Locally Linear Embedding. Yet LPP is linear and more crucially is defined everywhere in ambient space rather than just on the training data points [86]. Suppose we have a collection of data points of n-dimensional real vectors drawn from an unknown probability distribution. In increasingly many cases of interest in machine learning and data mining, one is confronted with the situation where n is very large. However, there might be reason to suspect that the “intrinsic dimensionality” of the data is much lower. This leads one to consider methods of dimensionality reduction that allow one to represent the data in a lower dimensional space. LPP is a linear dimensionality reduction algorithm which builds a graph incorporating neighborhood information of the data set. Using the notion of the Laplacian of the graph, it computes a transformation matrix which maps the data points to a subspace. This linear transformation optimally preserves local neighborhood information in a certain sense. The representation map generated by the algorithm may be viewed as a linear discrete approximation to a continuous map that naturally arises from the geometry of the manifold [87]. The generic problem of linear dimensionality reduction is the following. Given a set 𝒙1,𝒙2,…,𝒙𝑚 𝑖𝑛 𝑹𝑛, find a transformation matrix A that maps these m points to a set of points 𝒚1,𝒚2,…,𝒚𝑚 𝑖𝑛 𝑹𝑙 (𝑙≪𝑛), such that 𝒚𝑖 “represents” 𝒙𝑖, where 𝒚𝑖=𝐴𝑇𝒙𝑖. This method is of particular applicability in the special case where 𝒙1,𝒙2,…,𝒙𝑚 ∈𝑀 and 𝑀 is a nonlinear manifold embedded in 𝑹𝑛. LPP is a linear approximation of the nonlinear Laplacian Eigenmap. The algorithmic procedure is formally stated below: 1. Constructing the adjacency graph: Let G denote a graph with m nodes. An edge is put between nodes i and j if 𝒙𝑖 and 𝒙𝑗 are “close”. There are two variations: a. 𝜖-neighborhoods. [Parameter 𝜖∈𝑹]. Nodes i and j are connected by an edge if ‖𝒙𝑖−𝒙𝑗‖2<𝜖 where the norm is the usual Euclidean norm in 𝑹𝑛. b. 𝑘 nearest neighbors. [Parameter 𝑘∈𝑵]. Nodes i and j are connected by an edge if i is among k nearest neighbors of j or j is among k nearest neighbors of i. 62 2. Choosing the weights: Here, as well, we have two variations for weighting the edges. W is a sparse symmetric 𝑚×𝑚 matrix with 𝑊𝑖𝑗 having the weight of the edge joining vertices i and j, and 0 if there is no such edge. a. Heat kernel. [Parameter 𝑡∈𝑹]. If nodes i and j are connected, put 𝑊𝑖𝑗=𝑒−‖𝒙𝑖−𝒙𝑗‖2 𝑡 (3.40) The justification for this choice of weights can be traced back to [87]. b. Simple-minded. [No parameter]. 𝑊𝑖𝑗=1 if and only if vertices i and j are connected by an edge. 3. Eigenmaps: Compute the eigenvectors and eigenvalues for the generalized eigenvector problem: 𝑋𝐿𝑋𝑇𝒂=λ𝑋𝐷𝑋𝑇𝒂 (3.41) where D is a diagonal matrix whose entries are column (or row, since W is symmetric) sums of W, 𝐷𝑖𝑖=𝛴𝑗𝑊𝑗𝑖.𝐿=𝐷−𝑊 is the Laplacian matrix. The 𝑖𝑡ℎ column of matrix X is 𝒙𝑖. Let the column vectors 𝒂0,…,𝒂𝑙−1 be the solution of equation (3.41), ordered according to their eigenvalues, λ0<⋯<λ𝑙−1. Thus, the embedding is as follows: 𝒙𝑖→𝒚𝑖=𝐴𝑇𝒙𝑖,𝐴=(𝒂0,𝒂1,…,𝒂𝑙−1) (3.42) where 𝒚𝑖 is a 𝑙-dimensional vector, and A is a 𝑛×𝑙 matrix. An experiment was conducted with the Multiple Features Database [88]. This dataset consists of features of handwritten numbers (‘0’-‘9’) extracted from a collection of Dutch utility maps. 200 patterns per class (for a total of 2,000 patterns) have been digitized in binary images. These data points are mapped to a 2-dimensional space using different dimensionality reduction algorithms, PCA, LPP, and Laplacian Eigenmaps. The experimental results are shown in Figure 3.17 [86]. 63 Figure 3.17: The handwritten digits (‘0’-‘9’) are mapped into a 2-dimensional space. Each color corresponds to a digit. LPPs are obtained by finding the optimal linear approximations to the eigenfunctions of the Laplace Beltrami operator on the manifold. As a result, LPP shares many of the data representation properties of nonlinear techniques such as Laplacian Eigenmap. However, LPP is computationally much more tractable. Moreover, LPP is derived by preserving local information; hence it is less sensitive to outliers than PCA and it can has more discriminating power than PCA. 3.1.3.1.5. Neighborhood Preserving Embedding (NPE) Neighborhood Preserving Embedding (NPE) is a subspace learning algorithm different from Principal Component Analysis (PCA) which aims at preserving the global Euclidean structure, NPE aims at preserving the local neighborhood structure on the data manifold. Therefore, NPE is less sensitive to outliers than PCA. Also, NPE is defined everywhere, rather than only on the training data points. Furthermore, NPE may be 70 Figure 3.20: Output after kernel PCA with 𝑘(𝑥,𝑦)=(𝑥𝑇𝑦+1)2. The three groups are distinguishable using the first component only. Now consider a Gaussian kernel: 𝑘(𝑥,𝑦)=𝑒−‖𝑥−𝑦‖2 2𝜎2 (3.62) That is, this kernel is a measure of closeness, equal to 1 when the points coincide and equal to 0 at infinity. 71 Figure 3.21: Output after kernel PCA, with a Gaussian kernel. Note in particular that the first principal component is enough to distinguish the three different groups, which would be impossible using only linear PCA, because linear PCA operates only in the given (in this case two-dimensional) space, in which these concentric point clouds are not linearly separable. 3.1.3.2. Spectral unmixing Techniques In many studies, hyperspectral analysis techniques are divided into full-pixel and mixed-pixel classification techniques [96][97][98], where each pixel vector defines a spectral signature or fingerprint that uniquely characterizes the underlying materials at each site in a scene. Full-pixel classification techniques assume that each spectral signature comprises the response of one single underlying material. Often, this is not a realistic assumption, because the spatial resolution of the sensor is not fine enough to separate different pure signature classes at a macroscopic level. In consequence, these can jointly occupy a single pixel, and the resulting spectral signature will be a composite of the 72 individual pure spectral, often called endmembers in hyperspectral imaging terminology [99]. In this chapter, we explore the use of spectral unmixing for feature extraction prior to supervised classification of hyperspectral data using SVM. Moreover, differently from most feature extraction techniques available in literature, the features obtained using linear spectral unmixing are potentially easier to interpret due to their physical meaning [100]. Hyperspectral data are often used to determine what materials are present in a scene. In our case, materials of interest could include some tissues as tumor, healthy, necrosis in addition to blood, veins and surgical elements, etc. In fact, each pixel of a hyperspectral image could be compared to a material database to determine the type of material making up the pixel. Spectral unmixing aims at the decomposition of the mixed pixel spectrum into its constituent spectra, also called endmembers [101]. Each pixel in the hyperspectral image can be considered as being composed of linear combination of ground spectra or endmembers with each endmember contributing to the pixel spectra. Thus the spectral signature at each pixel in a L-dimensional hyperspectral image, 𝑌∈𝑅𝐿, when p is the number of endmembers, can be expressed as, 𝑦=𝑥+ 𝑛 (3.63) where, y - L-dimensional pixel vector. x and n - L-dimensional signal and noise vectors respectively. Since the signal vectors lie in an unknown p-dimensional subspace, each signal vector is given as, 𝑥=𝑀𝑠=∑𝑚𝑖𝑠𝑖 𝑝 𝑖=1 (3.64) where, 𝑀 is a L×p matrix, whose columns are L×1 endmembers. 𝑠 is the abundance fraction of each endmember in a pixel. 73 In essence, spectral unmixing can be defined as the process of determination of the number of image endmembers and their pure signatures and the amount in which they appear in the given mixes pixel [58]. Namely, Spectral unmixing consists of estimating the fraction of the pixel area covered by each material present in the scene [4]. Pixel unmixing algorithms can be separate into two main areas: Endmember Determination and Abundance Estimation algorithms as shown in figure 3.22. Figure 3.22: Spectral unmixing Diagram Process. Before describing the algorithms considered, it should emphasize the context in which these algorithms will be applied. In this regard, figure 3.23 [103] describes the classic method for analyzing hyperspectral images using linear mixed model. As shown in figure 3.23 [103], the methodology starts from a hyperspectral image and performs the following steps: 1. Dimensionality reduction: This step is optionally used by certain algorithms to reduce the computational load of successive steps by removing noise and redundant information in the image. 2. Endmembers extraction: In this step the pure spectral signatures that combine to result in mixed pixel in the image are identified. 3. Abundances estimation: The abundance of pure spectral signatures or endmembers is estimated at each pixel of the image. 74 Figure 3.23: classic method for analyzing hyperspectral images using linear mixed model. The spectral unmixing chain can be divided in three steps: number of endmember estimation, endmember extraction and abundances estimation. 3.1.3.2.1. Number of endmember estimation First, the number of endmembers 𝑝 is estimated directly from the original ndimensional hyperspectral image. For this purpose, two standard techniques widely used in the literature such as the HySime method [59] and the VD concept [104] are recommended. These techniques are described in the next section 3.1.4 corresponding to feature selection. In our case, the estimation of number of endmembers is not necessary because this number is known to us a priori and associated with the number of classes. 75 3.1.3.2.2. Endmember extraction algorithms Once the number of endmembers 𝑝 has been estimated, an automatic algorithm to extract a set of endmembers from the original hyperspectral image is applied [105]. There are many algorithms to unmix hyperspectral data each with their own strengths and weaknesses. Many algorithms assume that pure pixels (pixels which contain only one material) are present in a scene. Endmember extraction algorithms can be divided into two groups: first those based only in image spectral information and secondly those that uses both spectral and spatial information. In this dissertation we have considered only the algorithms in the first group, where the used endmember extraction techniques have been: OSP [106], VCA [107] and N-FINDR [108] as representative automatic algorithms of the trend based on use only spectral information. Some examples of algorithms belonging to the second group are: Automatic Morphological Endmember Extraction (AMEE) [109], Spatial Spectral Endmember Extraction (SSEE) [110] y Spatial Pre-Processing (SPP) [111] as representative automatic algorithms of the trend based using both spatial and spectral information. The endmember extraction algorithms based on spectral information depends on the ability to discriminate the pixels based solely on its spectral features. The following sections describe, the OSP, N-FINDR and VCA algorithms as three of the most representative techniques within this category. Other classic techniques such as PPI algorithm have been discarded for this study because of its semi-supervised nature, making it difficult to obtain consistent results in several automated executions. 3.1.3.2.2.1. Orthogonal Subspace Projection (OSP) The OSP algorithm was initially developed to find spectral signatures using the concept of orthogonal projections [112][113]. The algorithm uses an orthogonal projection operator that is given by the expression [114]: 𝑃𝑈⊥=𝐼−𝑈(𝑈𝑇𝑈)−1𝑈𝑇 (3.65) where 𝑈 is a matrix of spectral signatures, 𝑈𝑇 is the transpose of this matrix, and 𝐼 is the identity matrix. The algorithm uses the operator shown in equation (3.65) repetitively until finding a set of p orthogonal pixels from an initial pixel. The iterative process carried out by this algorithm can be summarized as follows: 76 1. Calculate 𝑡0, the brightest pixel in the hyperspectral imaging, using the following expression: 𝑡0=𝑎𝑟𝑔{𝑚𝑎𝑥(𝑥,𝑦)𝐹(𝑥,𝑦)𝑇∙𝐹(𝑥,𝑦)}, where 𝐹(𝑥,𝑦) is the pixel in the (𝑥,𝑦) coordinates of the image. As can be seen, the brightest pixel is the one with the greatest value when performing the cross product between the vector associated with that pixel and its transposed 𝐹(𝑥,𝑦)𝑇, or what is the same, the first norm of the pixel. 2. Applying an orthogonal projection operator denoted as 𝑃𝑈⊥, based on the expression (3.65), with 𝑈=𝑡0. This operator is applied to all pixels in the hyperspectral imaging. 3. Then the algorithm finds a new endmember called 𝑡0 with the highest value in the complementary space 〈𝑡0〉⊥, orthogonal to 𝑡0, as follows: 𝑡1=𝑎𝑟𝑔{𝑚𝑎𝑥(𝑥,𝑦)[𝑃𝑈⊥∙𝐹(𝑥,𝑦)]𝑇∙[𝑃𝑈⊥∙𝐹(𝑥,𝑦)]}. In other words, the algorithm searches for the pixel with the highest orthogonality with respect to 𝑡0. 4. The next step is to modify the 𝑈 matrix by adding the new endmember found, i.e., 𝑈=[𝑡0𝑡1]. 5. Next, the algorithm finds a new endmember called 𝑡2 with the highest value in the complementary space 〈𝑡0,𝑡1〉⊥, orthogonal to 𝑡0 and 𝑡1, as follows: 𝑡2=𝑎𝑟𝑔{𝑚𝑎𝑥(𝑥,𝑦)[𝑃𝑈⊥∙𝐹(𝑥,𝑦)]𝑇∙[𝑃𝑈⊥∙𝐹(𝑥,𝑦)]}. It should be noted that, unlike step 3) in which 𝑈=𝑡0, at this point the orthogonal projector is based on a matrix 𝑈=[𝑡0𝑡1]. 6. The process is repeated iteratively until the desired number of endmembers is found. This algorithm is effective in identifying a set of spectrally differentiated endmembers by the orthogonality condition imposed in the search process. As a negative feature, the algorithm can be sensitive to outliers and anomalous pixels [103][106]. 3.1.3.2.2.2. Vertex Component Analysis (VCA) This algorithm also makes use of the concept of orthogonal projections. However, unlike the previously described algorithm OSP, VCA algorithm exploits the fact that the endmember are the vertices of a simplex and that the affine transformation of a simplex is also a simplex. As a result, VCA models the data using a positive cone whose projection in an appropriately chosen hyperplane is another simplex whose vertices are the final 77 endmembers. After projecting the data in the selected hyperplane, the VCA algorithm projects all the image pixels in a random direction and uses the pixel with the greatest projection as the initial endmember. The others endmembers are sequentially identified by projecting iteratively the data in a direction orthogonal to the subspace spanned by the endmembers currently selected. The new endmember is then selected as the corresponding pixel to the most extreme projection, and the procedure is repeated until p endmembers are selected [103][107]. 3.1.3.2.2.3. N-FINDR N-FINDR algorithm [115][116] uses a technique based on identifying endmembers as the vertices of the simplex with largest volume that can form on the set of points. NFINDR does not work with the whole data hypercube but with a simplification of it with as many bands as endmembers are desired to find. For this type of dimensionality reduction are often used techniques as explained above, such as Principal Component Analysis (PCA) [117], Minimum Noise Fraction (MNF) [117] or Independent Component Analysis (ICA) [82]. The only parameter that has this algorithm is the number of endmembers to identify. The operation of the algorithm is described in the following steps: 1. Make a reduction of the image to a number of bands equal to the number of endmembers to be extracted, using one of the dimensionality reduction algorithms discussed above. 2. Select a random number of pixels that are labeled as endmembers. This initial selection will be refined iteratively. 3. The third step is to select a pixel of the original image. This pixel will successively exchange with each endmembers initially selected. 4. As the pixel will exchange with the initial endmembers, the volume of the hyperpolygon formed with the new point considered is calculated. 5. If the volume obtained after the exchange is greater than it had before the exchange, the new point results in a replacement in the endmembers set and the new pixel becomes part of endmembers set. Otherwise, the exchange is reversed. 6. Steps 3-5 are iteratively repeated to check all pixels of the image. Of such form that at the end of process we have a set of endmembers such that its volume is the largest possible. 78 It should be noted that in the second step of the algorithm, an initial set of endmembers is randomly established. If the initial estimate is suitable, the algorithm will reach the optimal solution. On the contrary, an erroneous initial estimate may result in failure to reach the optimal solution but we stay in a local maximum of the growth function of the hyperpolygon [116]. The algorithm assumes that an increase in the volume of the defined hyperpolygon by incorporating a new pixel in the endmembers set carries a higher quality thereof. However, figure 3.24 shows that the fact of using a polygon of greater volume does not ensure a better description of the set of points [103]. A more reliable parameter is the increase in the number of pixels that can be described using the new set of endmembers. Figure 3.24: N-FINDR algorithm operation. To conclude the description of this method, it is important that the endmembers identified by the algorithm N-FINDR correspond to pixels belonging to the original dataset. Using this algorithm, it is not possible to generate artificial endmembers because the replacements are always performed using existing points in the set of available samples. Thus, it may happen that the selected endmembers are not the most pure. Furthermore, the method is sensitive to outliers [103][108]. 79 3.1.3.2.3. Abundance estimation algorithms Finally, linear spectral unmixing (either unconstrained or constrained) can be used to estimate the abundance of each endmember in each pixel of the scene, providing a set of p abundance maps [100]. 3.1.3.2.3.1. Linear Spectral Unmixing In hyperspectral images, spectral mixing is the result of mixing of two or more spectrally distinct substances. Spectral unmixing is the process by which we can identify the constituents of the mixed pixel and their proportions. Generally, two models of mixing are assumed: linear and nonlinear [58]. As opposed to nonlinear unmixing, which generally require detailed information about physical properties that may not be always available, linear spectral unmixing consists of identifying the pure spectral components or endmembers. When the pure spectral signatures are identified, the proportion of each material in each pixel can be estimated. Abundances provide additional information about the composition of each pixel; if this information is used in a correct way, it may complement the results provided by traditional “hard” classification techniques. Moreover, non-linear mixing results from multiple scattering often due to non-flat surface [100]. The simplest and the most commonly assumed model for a mixed spectrum is a linear model [96]. A single pixel can be portrayed as a checkerboard mixture, assuming that there is no multiple scattering between components, then the spectral response of the pixel is a linear combination of the fractional abundances (area covered by each endmember in the pixel) of the individual substances [101], hence the term Linear Mixture Model (LMM). If there are 𝑝 endmembers, then the linear mixture model can be expressed as 𝑥=∑𝑚𝑖𝑠𝑖𝑗+𝑤𝑗=𝑀𝑠+𝑤 𝑝 𝑖=1 (3.66) 𝑗=1,2,……,𝑁 where, 𝑥 is the 𝐿×1 received pixel spectra. 𝑀 is a 𝐿×𝑝 matrix, whose columns are L×1 endmembers. 86 his work team at AT&T Labs [133]. It combines following ideas: the optimal hyperplanes search technique as a solution, the idea of convolution scalar product, the extension of the linear functions to nonlinear and the notion of soft margin to allow for errors in the training patterns. An important advantage is that it works on the principle of Structural Risk Minimization (SRM), it is better than ERM which use others many techniques. SVM then enables better generalization rather than a better classification of the training set (at the level of errors). There are two other reasons that have increased the interest in this newfangled classifier. SVM can be reduced to a problem of convex quadratic programming (QP), which is easier to solve compared to traditional methods and that seem to have a better performance (more robust) with high amounts of data. SVM has been used hitherto in many fields, such as: Text categorization, recognition of hand-written texts, image classification, bioinformatics, remote sensing, and now also in medical hyperspectral imaging, where it seems to have a higher performance than other classical techniques used [134]. 3.2.1.1 Theoretical foundations of SVM classifier This classifier belongs to the family of linear classifiers that induce linear separators or hyperplanes in spaces of high dimensionality characteristics, though may be easily adapted to act as non-linear classifiers by applying a function or not linear kernel in the input data. Its main objective is to get one surface (or hyperplane) capable of separating the different classes that can be grouped in a data distribution of a N-dimensional space, using an optimization process based on obtaining vectors which defining the class boundaries. These vectors are usually referred support vectors [135]. If we see the input data as two sets of vectors in a N-dimensional space, the SVM algorithm objective is simply to build a separating hyperplane in that space, which maximizes the margin of distance to the two data sets [136]. 87 Figure 3.27: Functional diagram of the SVM classifier. In figure 3.27 [5] we can see how to calculate this separating hyperplane, constructing two parallel hyperplanes, one on each side of the first. The two parallel hyperplanes are pushed or widened to as close as possible to the datasets. Intuitively, a good separation is achieved when the separating hyperplane is farthest from both classes. The greater the separation distance, the lower will generally be the classifier error [135]. In mathematical terms, given a training set of the equation: 𝐷={(𝑥𝑖,𝑐𝑖)|𝑥𝑖∈𝑅𝑝,𝑐𝑖∈{−1,1}}𝑖=𝑛 𝑛 (3.69) Where 𝑐𝑖 is 1 or -1, indicating the class to which 𝑥𝑖 belongs. Each 𝑥𝑖 is a real vector p-dimensional, we want to obtain a maximum distance hyperplane to the training sets which divide those belonging to 𝑐𝑖=1 of those with the value 𝑐𝑖=−1. Any hyperplane can be written as a set of points 𝑥 that satisfies the equation: 𝑤∙ 𝑥− 𝑏=0 (3.70) The 𝑤 vector is a normal vector perpendicular to the hyperplane. The 𝑏 ‖𝑤‖ parameter determines the displacement of the hyperplane over the origin. We want to choose the 𝑤 and the 𝑏 which maximize the distance between the two parallel hyperplanes, which are as far apart as possible depending on the data. These hyperplanes can be described with the formulas described below: 𝑤∙ 𝑥− 𝑏=1 (3.71) 88 and, 𝑤∙ 𝑥− 𝑏=−1 (3.72) Note that if the training set is linearly separable, we can choose two hyperplanes on the edge of the sets such that there are no points between them and then try to maximize their distance. Using geometry, we find that the distance between them is 2 ‖𝑤‖, so it is intended to minimize ‖𝑤‖. As we have to prevent the points are located in the boundary area, we add the constraint of equation 𝑤∙𝑥𝑖− 𝑏≥1 (3.73) for 𝑥𝑖 belonging to the first class and the restriction of equation 𝑤∙𝑥𝑖− 𝑏≤−1 (3.74) for 𝑥𝑖 in the second class. This can be written as: 𝑐𝑖(𝑤∙𝑥𝑖−𝑏)≥1,∀ 1≤𝑖≤𝑛 (3.75) We can compact the expression to reach the optimization problem: Choose 𝑤, 𝑏 to minimize ‖𝑤‖: Subject to 𝑐𝑖(𝑤∙𝑥𝑖−𝑏)≥1,∀ 1≤𝑖≤𝑛 (3.76) The optimization problem presented above is difficult because only depends on a value |𝑤|. The reason is that it is a non-convex optimization problem, which is known to be much more difficult to solve than the convex optimization problem. Fortunately, is possible to replace ‖𝑤‖ by 1 2‖𝑤‖2 without changing the solution. This is an optimization problem of quadratic programming. More clearly, the optimization problem can be reformulated as follows: Minimize 1 2‖𝑤‖2, subject to 𝑐𝑖(𝑤∙𝑥𝑖−𝑏)≥1,∀ 1≤𝑖≤𝑛 (3.77) The 1/2 factor is used as a mathematical convenience. Now, the problem presented to us can be resolved through programs and standard quadratic programming techniques. Writing the classification rule dual in its extended form, this reveals that the maximum distance to the hyperplane, and therefore the task of classification, it is only a 89 function of support vectors, ie, data that are on the borderline. The second form of SVM can be derived as the following expression: 𝑚𝑎𝑥∑𝛼𝑖 𝑛 𝑖=1 −12∑𝛼𝑖𝛼𝑗 𝑖,𝑗 𝑐𝑖𝑐𝑗𝑥𝑖𝑇𝑥𝑗 (3.78) Subject to 𝛼𝑖≥0, and ∑𝛼𝑖𝑐𝑖=0 𝑛 𝑖=1 (3.79) Where α terms are another representation of the weight vector in terms of the training set: 𝑤=∑𝛼𝑖𝑐𝑖 𝑖𝑥𝑖 (3.80) The original algorithm specifies a linear classifier, however, it may be modified to solve nonlinear classification problems, replacing the scalar product by a nonlinear kernel function. This allows the algorithm to determine the maximum distance to the hyperplane in a transformed features space. The transformation could be a nonlinear transformation and the transformed space could be a high dimensional space; this way, even though the classifier is a hyperplane in the space of high-dimensional features, it may not be linear in the original input space [137]. If the kernel used is Gaussian radial basis function type, the characteristics space is a Hilbert space of infinite dimension. The maximum distance is regulated, so the infinite dimension will not spoil the results [138]. Some of the commonly used kernels in SVM classifiers are listed below:  Polynomial (homogeneous): 𝑘(𝑥,𝑥′)=(𝑥∙𝑥′)𝑑  Polynomial (heterogeneous): 𝑘(𝑥,𝑥′)=(𝑥∙𝑥′+1)𝑑  Radial Basis Function: 𝑘(𝑥,𝑥′)=𝑒𝑥𝑝(−𝛾‖𝑥−𝑥′‖2),𝑓𝑜𝑟 𝛾>0  Gaussian Radial Basis Function: 𝑘(𝑥,𝑥′)=𝑒𝑥𝑝(‖𝑥−𝑥′‖2 2𝜎2)  Sigmoide: 𝑘(𝑥,𝑥′)=𝑡𝑎𝑛(𝑘𝑥∙𝑥′+𝑐),𝑓𝑜𝑟 𝑠𝑜𝑚𝑒 𝑘>0 𝑎𝑛𝑑 𝑐<0 In literature, we can find examples of kernels based on spectral metrics commonly used in hyperspectral analysis [139]. As stated at the beginning, processing used by the SVM does not require a large number of training patterns, as long as the chosen patterns are truly representative [140]. 90 3.2.1.2 Kernel functions available Our classifier has the following kernel functions that allow us to project information to a higher dimensional space which increases the computational capacity of the classifier:  Linear: 𝑢′∗ 𝑣 (3.81) Figure 3.28: The decision boundary of a Linear SVM.  Polynomial: (𝛾∗𝑢′∗𝑣+𝑐𝑜𝑒𝑓0)𝑛 (3.82) Figure 3.29: The decision boundary with a Polynomial kernel. 91  Radial basis function: 𝑒(−𝛾∗|𝑢−𝑣|2) (3.83) Figure 3.30: The decision boundary with a Radial Basis Function (RBF) kernel.  Sigmoid: tanh(𝛾∗𝑢′∗𝑣+𝑐𝑜𝑒𝑓0) (3.84) Figure 3.31: The decision boundary with a Sigmoid kernel.  Precomputed kernel: Allows us to introduce our own kernel using a matrix. By default, the classifier works with the radial basis function kernel. 92 3.3. Summary To take advantage of the large amount of information offered by hyperspectral imaging is necessary to design a robust preprocessing chain. This is because the high dimensionality that provides hyperspectral data leads to the presence of noise and redundant information which complicate its direct processing. To address this issue, preprocessing chains usually consist of two main parts: feature extraction and feature selection, in which it is possible to apply many different techniques. Firstly, feature extraction is responsible for extracting the main features of the hyperspectral data using transformations theory. Then feature selection is responsible of choosing the most important characteristics of this new set of transformed data. The result is an alternative representation of the original data set but with lower dimensionality and therefore much easier to process. The ultimate goal of this process is to classify the data in the most efficient way possible. For this purpose a Support Vector Machine (SVM) will be used. SVM is a very robust classifier that provides very good results working with this type of data. 93 Chapter 4 Proposed process for cancer detection using hyperspectral images In recent years the hyperspectral sensors technology has rapidly advanced providing data with very high spectral fidelity compared to multispectral systems. While these sensors facilitate the identification and classification, the high dimensionality and volume of data increases the bandwidth and computational complexity of the analysis. In addition, we have a large amount of redundancy in the hyperspectral data due to the high correlation between adjacent bands. In order to make an optimum classification and minimize the computational time, it is necessary to find a method to reduce the dimensionality of data while keeping the necessary spectral characteristics to classify data. The most general processing chain can be described with a small group of elaborate block chains. In fact, after obtaining a sample set (forming the hypercube, filtering the data, extracting samples from its contour, labeling samples, etc.) the relevant information can be extracted from the hyperspectral data. Starting from the sample set, three main steps can be described in the processing chain: transformation and feature extraction, feature selection and classification. 94 In this dissertation, feature extraction is implemented in two different ways to reduce the dimensions of the data: dimensionality reduction and unmixing techniques. The feature selection process enables efficient and selective process for choosing the best subset, based on the high-dimensional data. The combination of these features is complex, so fast and efficient algorithms are needed for brain cancer detection. Finally, the classification step produces a high level result ready for user interpretation. In this regard, each process contains a classification stage which is performed by the supervised pattern recognition algorithm called Support Vector Machine (SVM). In the beginning we only had ex vivo data available and in consequence the whole processing chain was used with this dataset while we waited to get the in vivo data. In this regard, we tested extensively all the processing chains with ex vivo data in order to select the processing chains that provided best results. Thus, when we had in vivo data available we could test the best combinations of processing with these new data. In this chapter, we present all the possible combinations of block chains selected for experimentation. The reason for the choice of the processing block chains that we will present below is based on previous studies in other research frameworks (some outside the medical field) with hyperspectral imaging, in which the techniques presented in this dissertation obtained good results. These chains showed good performance in terms of accuracy percentage in classification with respect to others. However, there is a world of possibilities where is possible to replace, add, or remove a block chain, or research on new techniques in order to obtain good results in characterizing the image. Based on previous experiments with different types of images or datasets and using different percentages of training sets, we have considered the chains that are described below. 95 4.1. Process #1: Support vector machine classifier used in brain cancer detection In this process the original data set is classified using SVM classifiers. In order to obtain the best possible results, the parameters of SVM classifier are studied. For this, all kernels available for this classifier are used to check the experimental results. The optimal parameters were selected using 10-fold cross-validation (selected after testing different configurations). Furthermore, it is shown the process of obtaining the set of samples from the files generated by the camera. Along this processing, the hyperspectral data cubes are formed from the files produced by the camera. These cubes have areas with samples of interest (healthy, tumor and necrosis tissues) and other areas (table, test plates and glow in the samples) which must be removed from the image. The hypercubes are spatially reduced in order to process it more easily and quickly, and for extracting samples removing the other elements in the image a contour functions are used. In addition, the original data are noisy due to the circumstances in which they were taken. For this reason, in this chain is intended to apply the types of filtering discussed above, in order to further improve the results of the classification. Finally all samples obtained from each hypercube are united in a single set of samples whose samples are properly labeled. Once obtained the sample set, it proceeds to the classification of the pixels using the SVM classification algorithm from the set of samples obtained and the training set. This stage produces a classified image with each pixel assigned to a particular class. This last step is common in almost all chains, so in the following sections where other processing chains are explained, it is possible to obviate the description of this block. In conclusion, taking as a decision criterion the best results obtained, this chain is intended to get two objectives: to establish the best possible configuration of SVM algorithm that will be used throughout this dissertation as well as to establish, depending on the type of filter used, the set of samples that will be used in other processing chains, while showing the process of obtaining the samples to classify. For this reason, the rest of processing chains begin directly from the set of samples obtained. The diagram used in this processing chain is shown in figure 4.1: 102 Figure 4.5: General flowchart of the processing chain #5. 4.6. Process #6: Patient simulation using SVM in brain cancer detection This process is based on dividing into two subsets of samples the ex vivo dataset as a final step before starting to work with in vivo samples, in order to simulate a real surgery. This is because in surgeries it is proposed to work on a database generated from generic spectral signatures obtained from different patients, so that a patient is classified using as training other patients from the database. Therefore by dividing the ex vivo dataset in two subsets is intended to simulate this situation, a subset of samples is classified obtaining the training set from the other subset of samples. The SVM classifier works extracting the training set from the data set that is to be classified. The SVM algorithm must be modified to work classifying a dataset and extracting the training of a different dataset. This modification is carried out in this chain to perform these tests and it will also be required to work with in vivo samples in which a patient is classified and the training is obtained from other patients. The subsets are totally independent, i.e., once separated, these are not joined at any time and these are divided into different sizes to perform several tests. Once again, in this processing chain the data sets from previous processing chains which obtained the best results are used. The diagram used in this processing chain is shown below in figure 4.6: SVM Random-Guided training Test samples Original, reduced or abundances set Test classification accuracy Random training Guided training 103 Figure 4.6: General flowchart of the processing chain #6. 4.7. Process #7: Consolidation of processing chains tested using in vivo data in brain cancer detection In this latest proposed process is intended to test the best combinations of techniques used in previous processing chains for ex vivo data now using in vivo data, checking how the change of data will affect the results. For working on this chain we have 5 different patients listed as 4, 5, 7, 8 and 10 referring to the number of surgical intervention. First, each patient is individually classified using the best combinations of techniques obtained in the previous processing chains. Test samples SVM Training samples Test classification accuracy Randomly selected Original, reduced or abundances set Classification set Training set 104 After that, we perform various combinations of patients, so we have patients to classify and patients to extract the trainings. In this regard several different tests will be performed in order to obtain the best procedure when establishing the training set facing real situations and always using the best processing chains obtained so far. In the case in which we use 2 groups of patients, one for classifying and one for training, these two datasets are independent and never are joined. For this reason, it is not possible to apply dimensionality reduction, because the transformation depends on the data set which is applied, then to apply the same technique of dimensionality reduction to two independent datasets generates that the data are transformed to different domains and these cannot be compared among them. For applying dimensional reduction techniques to two independent data sets would be necessary to transform one of the data sets and use the eigenvectors and eigenvalues obtained for transforming the second data set, so that both sets of data are transformed to the same domain. The diagram used in this processing chain is shown in figure 4.7: Figure 4.7: General flowchart of the processing chain #7. Test samples SVM Training samples Test classification accuracy Randomly selected Classifying patients Training patients 105 4.8. Summary In this section the different process developed in this dissertation are presented. All the processing chains are exhaustively tested using ex vivo samples. In summary, based on the results we should make the following decisions: 1. To decide whether to filter the data using the proposed filters or not to filter the data. 2. To decide whether to process the data using the dimensionality reduction techniques, unmixing techniques or a combination of both, or not to process the data. 3. To decide how to be extracted the training sets: randomly, guided or a combination of both. Once the best options are set, the resulting chains are tested with in vivo samples. 106 Chapter 5 Experimental results In this chapter we conducted an experimental evaluation of the different process of hyperspectral data described in the previous sections. Processing chains are evaluated according to different criteria, such as accuracy in the classification results obtained or its computational time. The main difficulty associated with thematic classification techniques is that usually there are several possible classes associated with different targets. The goal is ultimately to determine the presence or absence of each of the targets considered in each pixel, situation that can be expressed as a binary classification problem which can subsequently be extended to any number of classes. Several techniques are developed in order to evaluate the classification process in hyperspectral images. The general process is illustrated by a simple diagram in figure 5.1. As shown in this figure, the training process consists in using a subset of this information to train the supervised classifier (training patterns) and evaluate performance of the classifier with the rest of patterns labeled (test patterns), as it is shown in figure 5.1. 107 Figure 5.1: General process of a supervised system. Although in this project are used guided trainings, the distinction between training patterns and test patterns is generally performed randomly, trying to minimize the maximum number of training patterns necessary to achieve a satisfactory result of classification with the rest of samples, mainly because it is often difficult to obtain training patterns in this kind of applications. Due to the high dimensionality of the original data is attempted to compensate the need for a large number of samples using (optionally) data dimensionality reduction techniques. Moreover, due to the presence of mixed pixels in the hyperspectral imaging is attempted to solve the problem by applying (optionally) unmixing techniques. Therefore, throughout this chapter is intended to validate the combination of a set of techniques of dimensionality reduction and unmixing with a supervised classifier SVM. As for the sets of hyperspectral imaging selected to perform the experiments, these have been fully described in the corresponding section 2.2.3. However, it is worth noting that in this dissertation we have used two sets of samples: ex vivo and in vivo. Ex vivo samples are extracted from a real surgery after being removed from the patient, these are used to test all possible processing chains. When the best possible processing chains are established, in vivo samples obtained directly during the surgeries are used to check the performance in real situations of selected processing chains. Supervised classifier Training patterns Test patterns Hyperspectral imaging Randomly selected Test classification accuracy 108 5.1. Results of processing chains In this section we describe the results obtained by applying the different process presented in the previous chapter on the different hyperspectral imaging considered (ex vivo and in vivo), using ground-truth (this term refers to the accuracy of the training set's classification for supervised learning techniques) information to statistically validate the classification results obtained by applying the analysis methodology and evaluation metrics described above. Once the results obtained are presented separately for each processing chain, we proceed to discuss the results globally previously to the end of this chapter. It should be recalled that in all processing chains, except in which spectral unmixing techniques are applied, a supervised classifier SVM is used, so that the information that varies from one processing chain to another is mainly based on the preprocessing stage used. The classification is done by assigning a specific class label to each pixel once the learning phase of SVM classifier is completed, in which involves a number of training patterns selected from labeled pixels in the image from ground-truth information. Followed we present a detailed quantitative and comparative study analyzing the results obtained by different processing chains of hyperspectral imaging used. 5.1.1. Results of support vector machine classifier used in brain cancer detection In this experiment, we use the ex vivo samples data sets to analyze the possible kernels available in the SVM and the need to filter the datasets to remove noise from the samples prior to classification. For this processing chain we generate ten training sets by randomly selecting 3%, 5% and 10% of the ground-truth pixels. Then, the three considered types of input samples (original, HySime filter and smooth filter) are built for the selected training samples and used to train an SVM classifier in which two types of kernels: linear and radial basis function (RBF) are used. The SVM was trained with each of these training subsets and then evaluated with the remaining test set. Each experiment was repeated ten times, and the mean and standard deviation of the 109 different evaluation metrics values were reported. Kernel parameters were automatically optimized and selected by the internal function of SVM. The rest of kernels offered by SVM were tested also, but linear and RBF kernels were selected for being the most representative. The followings tables summarizes the evaluation metrics obtained after applying the considered SVM classification system (with linear and radial basis function (RBF) kernels) to the set of samples obtained after applying the filters considered in the process #1 (see 4.1) to the ex vivo samples. To refer to the filters used and the classes present, the nomenclature used is the same as that set forth previously in the section 3.1.2. AA No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 88.08% (0.0101) 88.76% (0.0079) 89.35% (0.0065) 94.90% (0.0047) 95.61% (0.0034) 96.65% (0.0011) 96.57% (0.0037) 97.14% (0.0030) 98.27% (0.0018) RBF 77.27% (0.0259) 83.84% (0.0196) 90.78% (0.0119) 79.31% (0.0370) 87.65% (0.0177) 92.56% (0.0187) 82.28% (0.0268) 89.73% (0.0199) 95.71% (0.0127) OA No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 85.76% (0.0119) 86.59% (0.0082) 87.56% (0.0089) 93.53% (0.0066) 94.65% (0.0050) 95.91% (0.0014) 95.41% (0.0048) 96.23% (0.0039) 97.55% (0.0028) RBF 74.21% (0.0272) 80.90% (0.0212) 88.93% (0.0143) 77.10% (0.0405) 85.15% (0.0216) 91.08% (0.0199) 79.16% (0.0263) 87.75% (0.0220) 94.50% (0.0152) Table 5.1: Classification accuracies (percentage) and standard deviation obtained after applying the considered SVM classification system (with linear and RBF kernels) to three different types of input datasets (original, HySime filter and smooth filter) extracted from the ex vivo samples (ten randomly chosen training sets). 110 Sensitivity No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 92.02% (0.0110) 92.99% (0.0168) 94.72% (0.0086) 96.92% (0.0119) 97.30% (0.0077) 98.19% (0.0045) 98.75% (0.0044) 98.82% (0.0057) 99.73% (0.0015) RBF 75.95% (0.0585) 83.17% (0.0226) 90.80% (0.0238) 78.83% (0.0651) 87.21% (0.0351) 92.30% (0.0284) 81.25% (0.0404) 90.25% (0.0335) 96.53% (0.0116) Specificity No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 98.91% (0.0025) 99.07% (0.0021) 99.18% (0.0013) 99.64% (0.0011) 99.79% (0) 99.85% (0) 99.87% (0.0012) 99.98% (0) 99.98% (0) RBF 93.78% (0.0107) 96.40% (0.0113) 98.50% (0.0029) 94.43% (0.0183) 97.80% (0.0052) 98.85% (0.0045) 95.96% (0.0133) 98.37% (0.0050) 99.57% (0.0022) Table 5.2: Sensitivity and specificity for class 1 (percentage) and standard deviation obtained after applying the considered SVM classification system (with linear and RBF kernels) to three different types of input datasets (original, HySime filter and smooth filter) extracted from the ex vivo samples (ten randomly chosen training sets). Sensitivity No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 94.22% (0.0078) 94.67% (0.0074) 95.10% (0.0035) 97.97% (0.0021) 98.18% (0.0046) 98.69% (0.0011) 98.86% (0.0028) 99.08% (0.0029) 99.57% (0.0013) RBF 88.10% (0.0124) 92.51% (0.0123) 96.11% (0.0059) 88.95% (0.0219) 94.58% (0.0098) 97.00% (0.0100) 91.58% (0.0185) 95.56% (0.0093) 98.60% (0.0052) Specificity No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 77.35% (0.0169) 78.19% (0.0134) 78.95% (0.0145) 88.32% (0.0125) 90.27% (0.0139) 92.08% (0.0034) 91.03% (0.0094) 92.45% (0.0094) 94.59% (0.0070) RBF 64.66% (0.0287) 71.40% (0.0250) 81.04% (0.0214) 68.37% (0.0459) 76.59% (0.0283) 84.39% (0.0292) 69.45% (0.0288) 79.86% (0.0305) 89.27% (0.0255) Table 5.3: Sensitivity and specificity for class 2 (percentage) and standard deviation obtained after applying the considered SVM classification system (with linear and RBF kernels) to three different types of input datasets (original, HySime filter and smooth filter) extracted from the ex vivo samples (ten randomly chosen training sets). 111 Sensitivity No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 61.86% (0.0257) 62.60% (0.0259) 62.23% (0.0223) 78.37% (0.0216) 81.74% (0.0259) 84.38% (0.0093) 81.77% (0.0188) 84.47% (0.0191) 87.85% (0.0155) RBF 49.81% (0.0291) 56.86% (0.0304) 68.79% (0.0246) 53.69% (0.0407) 63.29% (0.0287) 73.78% (0.0372) 54.42% (0.0380) 66.63% (0.0371) 79.61% (0.0445) Specificity No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 96.64% (0.0048) 96.86% (0.0054) 97.05% (0.0026) 98.82% (0.0017) 98.85% (0.0033) 99.16% (0) 99.31% (0.0015) 99.37% (0.0020) 99.71% (0) RBF 94.07% (0.0113) 96.29% (0.0059) 98.03% (0.0036) 94.45% (0.0102) 97.17% (0.0048) 98.41% (0.0054) 95.84% (0.0082) 97.67% (0.0067) 99.20% (0.0030) Table 5.4: Sensitivity and specificity for class 3 (percentage) and standard deviation obtained after applying the considered SVM classification system (with linear and RBF kernels) to three different types of input datasets (original, HySime filter and smooth filter) extracted from the ex vivo samples (ten randomly chosen training sets). Kappa No filter Filter 1 Filter 2 Kernel TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% TR. 3% TR. 5% TR. 10% Linear 0.7623 (0.0185) 0.7742 (0.0127) 0.7858 (0.0141) 0.8893 (0.0109) 0.9073 (0.0081) 0.9273 (0.0024) 0.9204 (0.0080) 0.9339 (0.0066) 0.9558 (0.0049) RBF 0.5829 (0.0399) 0.6849 (0.0321) 0.8086 (0.0231) 0.6250 (0.0618) 0.7521 (0.0336) 0.8448 (0.0335) 0.6587 (0.0410) 0.7920 (0.0355) 0.9024 (0.0260) Table 5.5: Kappa coefficient and standard deviation obtained after applying the considered SVM classification system (with linear and RBF kernels) to three different types of input datasets (original, HySime filter and smooth filter) extracted from the ex vivo samples (ten randomly chosen training sets). As shown by the tables, the classification results are correlated with the training set size (the larger the training set, the better the classification results and the computing times). The good generalization ability exhibited by SVM is demonstrated by the classification results reported for the original spectral information, even with very limited training sets. These results can be improved in this chain by two different ways: kernels used and filtering samples. On the one hand, regarding the kernels used we show a comparison between linear kernel and RBF kernel, from which we can see that the linear kernel performs better than RBF for this particular application. In addition to the better 118 Figure 5.4: Computing time for the different dimensionality reduction techniques applied. Figure 5.5: AA, OA and Kappa for the different dimensionality reduction techniques applied. 0 50 100 150 200 250 300 Original PCA18 PCA21 ICA18 ICA21 LPP18 LPP21 NPE18 NPE21 Time (sec) 0 20 40 60 80 100 120 AA (%) OA (%) Kappa*100 119 Figure 5.6: Sensitivity of class 1 for the different dimensionality reduction techniques applied. Figure 5.7: Specificity of class 1 for the different dimensionality reduction techniques applied. 84 86 88 90 92 94 96 98 100 102 Sensitivity class 1 75 80 85 90 95 100 105 Specificity class 1 120 Figure 5.8: Sensitivity for the different dimensionality reduction techniques applied. Figure 5.9: Sensitivity for the different dimensionality reduction techniques applied. As shown in the above tables and graphs, the best results applying dimensionality reduction is obtained by MNF. Nevertheless, MNF technique was applied using the ENVI software. This is proprietary software and it cannot be included in this dissertation because the ULPGC has not purchase a license. Therefore, the ICA technique is the best option in terms of dimensionality reduction because this offers very competitive results and reduces the computational cost. 0 20 40 60 80 100 120 Sensitivity Class 1 Class 2 Class 3 0 20 40 60 80 100 120 Specificity Class 1 Class 2 Class 3 121 However, the original information provides the better accuracy results in comparison to the rest of the chain but the computational cost increases significantly. An alternative way of applying dimensionality reduction is a special version from smooth, which is called smooth reduction filter in this dissertation. This version was developed by a coworker in the department and allows us to reduce the number of bands of the original dataset similar to dimensionality reduction. In this regards, we apply the following process chain recommended by smooth reduction filter developer: 1. Denoising the samples using HySime filter. 2. Apply smooth reduction filter to reduce to 115 bands. 3. Implement a dimensionality reduction stage. 4. Classify the dataset. Once the dataset is obtained after applying the smooth reduction filter, the feature selection step is applied, obtaining the results shown in the table 5.13: Technique Virtual Dimensionality V.D. mean HySime Features 22 21 20 20 18 20.20 26 Table 5.13: Optimal number of features to samples filtered by smooth reduction. In this case, the values obtained from different techniques are similar, so we have chosen the average value obtained by virtual dimensionality technique and the value obtained by HySime, namely 20 and 26. In the table 5.14 the results obtained of applying smooth reduction filter prior to dimensionality reduction stage are shown. By way of comparison, in the firs column of results in this table, the results obtained of classifying directly the dataset without dimensionality reduction stage are shown. 122 Smooth reduction filter Smooth reduction 115 PCA 20 PCA 26 ICA 20 ICA 26 AA 95.39% (0.0033) 92.46% (0.0039) 93.50% (0.0036) 92.54% (0.0045) 93.55% (0.0043) OA 94.33% (0.0044) 90.74% (0.0046) 92.18% (0.0049) 90.92% (0.0043) 92.34% (0.0050) Kappa 0.8999 (0.0075) 0.8393 (0.0074) 0.8633 (0.0080) 0.8422 (0.0072) 0.8658 (0.0083) Time (sec) 57.7579 51.7009 51.3960 89.4440 63.0582 Table 5.14: Results of applying smooth reduction filter prior dimensionality reduction and classification. To summarize, the results obtained are reported in the following bar graphs, from figure 5.10 to figure 5.12: Figure 5.10: Average Accuracy for the application of smooth reduction filter prior dimensionality reduction stage. 90,5 91 91,5 92 92,5 93 93,5 94 94,5 95 95,5 96 Smooth reduction 115 PCA20 PCA26 ICA20 ICA26 Average Accuracy (%) 123 Figure 5.11: Overall Accuracy for the application of smooth reduction filter prior dimensionality reduction stage. Figure 5.12: Results obtained for the application of smooth reduction filter prior dimensionality reduction stage. As in the first results, in this case the best results obtained were without dimensionality reduction stage. In this regard, we proceeded to check the rest of available filter combinations trying to improve results. The combinations proposed were: 88 89 90 91 92 93 94 95 Smooth reduction 115 PCA20 PCA26 ICA20 ICA26 Overall Accuracy (%) 0 20 40 60 80 100 120 HySimesmooth115 PCA20 PCA26 ICA20 ICA26 Smooth reduction filter AA (%) OA (%) Kappa*100 Time (sec) 124  Combination 1: First apply HySime filter and then apply smooth reduction filter to reduce to 115 bands (as in the previous case).  Combination 2: Directly apply the smooth reduction filter to reduce to 115 bands.  Combination 3: First apply original smooth filter and then apply smooth reduction filter to reduce to 115 bands. The results obtained for these tests are shown in the followings tables and bar graphs: Filtering order Combination 1 Combination 2 Combination 3 AA 95.39% (0.0033) 90.73% (0.0035) 96.96% (0.0020) OA 94.33% (0.0044) 89.02% (0.0064) 96.08% (0.0022) Kappa 0.8999 (0.0075) 0.8097 (0.0097) 0.9298 (0.0039) Time (sec) 57.7579 105.1719 44.7789 Table 5.15: Results of applying the different filter combinations prior classification. Figure 5.13: Average Accuracy for the different filter combinations prior classification. 87 88 89 90 91 92 93 94 95 96 97 98 Combination 1 Combination 2 Combination 3 Average Accuracy (%) 125 Figure 5.14: Overall Accuracy for the different filter combinations prior classification. Figure 5.15: Results of applying the different filter combinations prior classification. Finally, in order to find the best way to reduce the number of bands, a comparison between the different techniques seen so far is performed. In this regard, the following processes are proposed, combining the best results obtained right now:  Process 1: The dataset is filtered using original smooth, then smooth reduction filter is used to reduce to 115 bands and finally ICA with 21 features is applied. 84 86 88 90 92 94 96 98 Combination 1 Combination 2 Combination 3 Overall Accuracy (%) 0 20 40 60 80 100 120 Combination 1 Combination 2 Combination 3 Filtering order AA (%) OA (%) Kappa*100 Time (sec) 126  Process 2: The dataset obtained from process #1 is dimensionally reduced using ICA with 21 features.  Process 3: The dataset is filtered using original smooth, then smooth reduction filter is used to reduce to 21 bands. The results obtained for this process are shown in the following tables and bar graphs: Comparative Process 1 Process 2 Process 3 AA 94.16% (0.0065) 94.78% (0.0054) 85.93% (0.0022) OA 92.49% (0.0087) 93.34% (0.0080) 81.87% (0.0037) Kappa 0.8683 (0.0145) 0.8826 (0.0134) 0.6988 (0.0049) Time (sec) 10.3263 9.5911 20.7463 Table 5.16: Results of applying different ways to reduce the number of bands. Figure 5.16: Average Accuracy for different ways to reduce the number of bands. 80 82 84 86 88 90 92 94 96 Process 1 Process 2 Process 3 Average Accuracy (%) 127 Figure 5.17: Overall Accuracy for different ways to reduce the number of bands. Figure 5.18: Comparative for different ways to reduce the number of bands. Theoretically, applying dimensionality reduction the classification results should be improved by removing the latest bands or features in the transformed dataset, which have redundant information and noise. However, as shown by the above tables, dimensionality reduction prior to classification cannot lead to improved classification results with regards to the process chain #1 and the original spectral information. This indicates that for this type of data, all the bands present useful information even in the case of the last bands, 76 78 80 82 84 86 88 90 92 94 96 Process 1 Process 2 Process 3 Overall Accuracy (%) 0 10 20 30 40 50 60 70 80 90 100 Process 1 Process 2 Process 3 Comparative AA (%) OA (%) Kappa*100 Time (sec) 134 5. Spectral unmixing chain 5: Dimensionality reduction is applied to training dataset using ICA to reduce 21 bands, then the endmembers are extracted from the reduced training dataset and the abundances are estimated. Abundance estimate FLSU LSU Endmembers extraction NFINDR NFINDR Class 1 1546/4469 34.5939% 1762/4469 39.4272% Class 2 3956/10881 36.3570% 4203/10881 38.6270% Class 3 2488/3230 77.0279% 2415/3230 74.7678% Total success 7990/18580 43.0032% 8380/18580 45.1023% Table 5.20: Endmembers extracted from training data set with dimensionality reduction using unmixing techniques. 6. Spectral unmixing chain 6: One endmember is extracted from each class of the dataset and then the abundances are estimated. Abundance estimate FLSU LSU Endmembers extraction NFINDR VCA OSP NFINDR VCA OSP Class 1 992/4469 22.1974% 2023/4469 45.2674% 4469/4469 100% 1781/4469 39.8523% 1903/4469 42.5822% 4465/4469 99.9105% Class 2 2904/10881 26.6887% 5862/10881 53.8737% 32/10881 0.2941% 643/10881 5.9094% 6692/10881 61.5017% 47/10881 0.4319% Class 3 777/3230 24.0557% 1421/3230 43.9938% 4/3230 0.1238% 892/3230 27.6161% 1464/3230 45.3251% 8/3230 0.2477% Total success 4673/18580 25.1507% 9306/18580 50.0861% 4505/18580 24.2465% 3316/18580 17.8471% 10059/18580 54.1389% 4520/18580 24.3272% Table 5.21: Each endmember extracted from each class using unmixing techniques. 135 7. Spectral unmixing chain 7: Dimensionality reduction is applied to dataset using ICA to reduce 21 bands, then one endmember is extracted from each class of the dataset and the abundances are estimated. Abundance estimate FLSU LSU Endmembers extraction NFINDR VCA OSP NFINDR VCA OSP Class 1 1789/4469 40.0313% 3408/4469 76.2587% 148/4469 3.3117% 1720/4469 38.4874% 4014/4469 89.8188% 141/4469 3.1551% Class 2 7945/10881 73.0172% 10254/10881 94.2377% 22/10881 0.2022% 9239/10881 84.9095% 9409/10881 86.4718% 22/10881 0.2022% Class 3 1145/3230 35.4489% 1835/3230 56.8111% 3223/3230 99.7833% 951/3230 29.4427% 2914/3230 90.2167% 3224/3230 99.8142% Total success 10879/18580 58.5522% 15497/18580 83.4069% 3393/18580 18.2616% 11910/18580 64.1012% 16337/18580 87.9279% 3387/18580 18.2293% Table 5.22: Each endmember extracted from each class with dimensionality reduction using unmixing techniques. In general, the use of spectral unmixing in this kind of samples doesn’t offer competitive results. This is because the data set is very diverse and a whole class cannot be represented using only one endmember. In order to try to solve this problem and improve the results obtained it is intended to extract the endmembers from an average of several pixels because in this way, using multiple pixels to form the endmembers, the classes will be better represented. In this case, the results are accompanied with its corresponding abundance maps of each class, displaying the fractional amount of that class present at each pixel. This allows us to visually verify the correct behavior of processes tested. In these maps the pixels distribution is as follows:  Class 1: It is comprised between the pixels 0 and 223.  Class 2: It is comprised between the pixels 224 and 767.  Class 3: It is comprised between the pixels 768 and 929. 136 Taking into account that for the above unmixing chains the best results are obtained from dimensionality reduction dataset, the processing for these remaining chains is always the same: dimensionality reduction is applied to dataset using ICA to reduce 21 bands, then the endmembers are extracted from the reduced dataset and the abundances are estimated. In these cases, they are differentiated by its endmembers extraction technique. The chain proposed is the following one: 8. Spectral unmixing chain 8: 6 endmembers are extracted from each class, and then the means of these endmembers are calculated for each class resulting in 3 endmembers (one per class). The abundances are estimate using these 3 endmembers formed by the means of extracted endmembers. Abundance estimate FLSU LSU Endmembers extraction NFINDR VCA OSP NFINDR VCA OSP Class 1 4192/4469 93.8017% 1849/4469 41.3739% 2881/4469 64.4663% 3827/4469 85.6344% 1829/4469 40.9264% 3730/4469 83.4639% Class 2 7861/10881 72.2452% 3948/10881 36.2834% 80/10881 0.7352% 7634/10881 70.1590% 3920/10881 36.0261% 108/10881 0.9926% Class 3 1263/3230 39.1022% 2879/3230 89.1331% 3150/3230 97.5232% 1877/3230 58.1115% 2882/3230 89.2260% 3080/3230 95.3560% Total success 13316/18580 71.6685% 8676/18580 46.6954% 6111/18580 32.8902% 13338/18580 71.7869% 8631/18580 46.4532% 6918/18580 37.2336% Table 5.23: 6 endmembers extracted from each class and the mean of these endmembers are calculated for each class, resulting in 3 endmembers using unmixing techniques. 137  Abundances map using NFINDR in spectral unmixing chain 8: Figure 5.24: Abundance map using NFINDR in spectral unmixing chain 8.  Abundances map using VCA in spectral unmixing chain 8: Figure 5.25: Abundance map using VCA in spectral unmixing chain 8. 138  Abundances map using VCA in spectral unmixing chain 8: Figure 5.26: Abundance map using VCA in spectral unmixing chain 8. 9. Spectral unmixing chain 9: 6 endmembers are extracted from each class and the abundances are estimate using the 18 endmembers obtained, then the estimated abundances for each class are added resulting in 3 abundances. Abundance estimate FLSU LSU Endmembers extraction NFINDR VCA OSP NFINDR VCA OSP Class 1 4036/4469 90.3110% 4253/4469 95.1667% 4203/4469 94.0479% 1931/4469 43.2088% 3500/4469 78.3173% 1894/4469 42.3808% Class 2 6871/10881 63.1468% 6881/10881 63.2387% 2893/10881 26.5876% 4999/10881 45.9425% 8437/10881 77.5388% 2884/10881 26.5049% Class 3 1913/3230 59.2260% 3095/3230 95.8204% 2709/3230 83.8700% 1815/3230 56.1920% 2684/3230 83.0960% 2817/3230 87.2136% Total success 12820/18580 68.9989% 14229/18580 76.5823% 9805/18580 52.7718% 8745/18580 47.0667% 14621/18580 78.6921% 7595/18580 40.8773% Table 5.24: 6 endmembers extracted from each class and abundances are added per class resulting in 3 abundances using unmixing techniques. 139  Abundances map using NFINDR in spectral unmixing chain 9: Figure 5.27: Abundance map using NFINDR in spectral unmixing chain 9.  Abundances map using VCA in spectral unmixing chain 9: Figure 5.28: Abundance map using VCA in spectral unmixing chain 9. 140  Abundances map using OSP in spectral unmixing chain 9: Figure 5.29: Abundance map using OSP in spectral unmixing chain 9. 10. Spectral unmixing chain 10: Endmembers extraction is not implemented. Instead, each endmember is calculated as the average of all pixels in its class and then the abundances estimation is performed. Abundance estimate FLSU LSU Datatype Dimensionality reduction Original dataset Dimensionality reduction Original dataset Class 1 4293/4469 96.0618% 2237/4469 50.0559% 4027/4469 90.1096% 2037/4469 45.5807% Class 2 9292/10881 85.3966% 4723/10881 43.4059% 9375/10881 86.1594% 6144/10881 56.4654% Class 3 2956/3230 91.5170% 1836/3230 56.8421% 2918/3230 90.3406% 1618/3230 50.0929% Total success 16541/18580 89.0258% 8796/18580 47.3412% 16320/18580 87.8364% 9799/18580 52.7395% Table 5.25: Each endmember is calculated as the average of all pixels in its class using unmixing techniques. 141  Abundances map using dataset with dimensionality reduction in spectral unmixing chain 10: Figure 5.30: Abundance map using dataset with dimensionality reduction in spectral unmixing chain 10.  Abundances map using original dataset in spectral unmixing chain 10: Figure 5.31: Abundance map using original dataset in spectral unmixing chain 10. 142 In general, this new approach to extract the endmembers based on obtaining these as an average of several pixels improves the results obtained, especially in the last unmixing chain with dimensionality reduction. This is because endmember extraction techniques are generally sensitive to outliers and anomalies. The objective of this procedure is to ensure that the pixels are as representative as possible of each class. Although the results have improved with the latest tests, the results are still worse than in the case of applying only dimensional reduction or in the case of not reducing the size of the data in any way and work on the original data. In the next section we will try to improve the results using combined techniques of unmixing and classification. 5.1.4. Results of Mixed techniques Unmixing-SVM used in brain cancer detection In this section we apply mixed techniques Unmixing-SVM in order to improve the results obtained in the last section in which only unmixing techniques were applied. For this purpose, the process chains proposed in this section are exactly the same as those used in the previous section, except that at the end a classification step is performed by the SVM classifier, which has as input data the estimated abundances. Processing chains are as follows: 143 1. Unmixing-SVM chain 1: The endmembers are extracted from the dataset and then the abundances are estimated. Abundance estimate FLSU LSU Endmembers extraction NFINDR VCA NFINDR VCA AA 44.49% (0.0055) 54.86% (0.0022) 52.42% (0.0022) 53.39% (0.0023) OA 28.60% (0.0072) 47.46% (0.0041) 46.89% (0.0106) 48.65% (0.0066) Kappa 0.0847 (0.0065) 0.2428 (0.0025) 0.2114 (0.0044) 0.2317 (0.0038) Time (sec) 36.9339 25.3817 15.0297 9.4290 Table 5.26: Endmembers extracted from data set using mixed unmixing-SVM techniques. 2. Unmixing-SVM chain 2: The endmembers are extracted from the normalized dataset and then the abundances are estimated. Abundance estimate FLSU LSU Endmembers extraction NFINDR VCA NFINDR VCA AA 42.83% (0.0126) 54.10% (0.0032) 52.42% (0.0022) 53.39% (0.0023) OA 24.89% (0.0067) 46.71% (0.0027) 46.89% (0.0106) 48.65% (0.0066) Kappa 0.0581 (0.0117) 0.2365 (0.0021) 0.2114 (0.0044) 0.2317 (0.0038) Time (sec) 33.8021 25.2906 15.1259 9.5866 Table 5.27: Endmembers extracted from normalized data set using mixed unmixing-SVM techniques.