Target detection of hyperspectral imagery on FPGAs

TARGET DETECTION OF HYPERSPECTRAL IMAGERY ON FPGAS DETECCIÓN DE OBJETIVOS EN IMÁGENES HIPERESPECTRALES SOBRE FPGAS Juan Martín Bárez Alonso GRADO EN INGENIERÍA DE COMPUTADORES FACULTAD DE INFORMÁTICA UNIVERSIDAD COMPLUTENSE DE MADRID Trabajo Fin de Grado Septiembre 2020 Directores: Carlos González Calvo José Manuel Mendías Cuadros Resumen en castellano En los últimos años ha ocurrido un resurgimiento de la carrera espacial motivado especialmente por empresas comerciales. Sus aeronaves son equipadas con una multitud de sensores, siendo uno de ellos las cámaras hiperespectrales. Este tipo de cámaras toma imágenes en cientos de bandas espectrales diferentes, con el objetivo de proporcionar información sobre el terreno. A causa del gran tamaño de las imágenes hiperespectrales, estas son enviadas a la Tierra para su procesado, con el consecuente coste de transmisión y almacenamiento. Preferentemente estas imágenes deberían procesarse o comprimirse in situ para enviar solo una fracción de los datos obtenidos. Dados el entorno espacial y las características este tipo de algoritmos, las FPGAs o ASICs se postulan como un sistema óptimo para su implementación. Este trabajo presenta una implementación sobre FPGAs del algoritmo Reed-Xiaoli de detección de anomalías para imágenes hiperespectrales. Para su implementación se ha realizado un análisis de las operaciones del algoritmo, centrada en una versión en punto otante y otra en aritmética de enteros, y de las repercusiones que tienen ciertas decisiones con la precisión que se alcanza. Demostrando de esta manera cómo algoritmos complejos con operaciones en punto otante pueden ser ejecutados en FPGAs al transformarlos para utilizar aritmética de enteros. Palabras clave Imágenes hiperespectrales, Algoritmo RX, Aritmética de punto otante, Aritmética de enteros, Hardware recongurable, VHDL. Abstract In recent years there has been a resurgence in the space race, motivated especially by commercial companies. Their aircrafts are equipped with a multitude of sensors, one of them being hyperspectral cameras. This type of camera takes images in hundreds of dierent spectral bands, with the aim of providing information of the ground. Because of the large size of hyperspectral images, they are sent to ground stations for processing, with the consequent cost of transmission and storage. Preferably these images should be processed or compressed on site and only a fraction of the data obtained should be sent. Given the spatial environment and the characteristics of these types of algorithms, FPGAs or ASICs are postulated as an optimal system for their implementation. This work presents an FPGA implementation of the Reed-Xiaoli algorithm of anomaly detection for hyperspectral images. For its implementation, an analysis of the operations of the algorithm has been made, centered in a oating point version and another one in integer arithmetic, and of the repercussions that certain decisions have with the precision that is reached. Thus, demonstrating how complex algorithms with oating point operations can be executed in FPGAs by transforming them to use integer arithmetic. Keywords Hyperspectral images, RX algorithm, oating point arithmetic, integer arithmetic, recongurable hardware, VHDL Table of Contents Índice i List of Figures iii List of Tables iv 1 Introduction 1 1.1 Motivation and objectives . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Relatedwork ............................... 2 1.2.1 Recongurable hardware . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Hyperspectral imagery . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3 Anomalydetection ........................ 8 1.2.4 Comparing oating and xed point . . . . . . . . . . . . . . . 11 1.3 Projectplan................................ 12 2 Software Model 13 2.0.1 Mean, deviation and covariancce matrix . . . . . . . . . . . . 14 2.0.2 Inverse............................... 15 2.0.3 Matrix multiplication . . . . . . . . . . . . . . . . . . . . . . . 17 2.1 Adapting the arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.1 DSPBlocks............................ 18 2.1.2 Determining multiplier size . . . . . . . . . . . . . . . . . . . . 19 2.1.3 Improving accuracy . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Implementation 20 3.1 General overview of the system . . . . . . . . . . . . . . . . . . . . . 21 3.2 Descriptionbymodule .......................... 21 3.2.1 Control .............................. 21 3.2.2 Inverter .............................. 23 3.2.3 Meansubtract........................... 26 3.2.4 Matrix multiplication . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.5 Coordinate sorter . . . . . . . . . . . . . . . . . . . . . . . . . 29 i 4 Results 30 4.1 Recongurable platform . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Hyperpectral image datasets . . . . . . . . . . . . . . . . . . . . . . . 30 4.3 Adequacy of approximation . . . . . . . . . . . . . . . . . . . . . . . 32 4.3.1 Floatingpoint........................... 32 4.3.2 Fixedpoint ............................ 32 4.3.3 Results for HYDICE . . . . . . . . . . . . . . . . . . . . . . . 33 4.3.4 ResultsforAVIRIS........................ 33 4.4 Computational eciency . . . . . . . . . . . . . . . . . . . . . . . . . 36 5 Conclusions 38 ii List of Figures 1.1 Generic FPGA hardware architecture . . . . . . . . . . . . . . . . . . 3 1.2 Eciency comparison between CPU, GPU and FPGA for image processing................................... 4 1.3 Overview of an heterogeneous FPGA . . . . . . . . . . . . . . . . . . 5 1.4 Comparison between a standard and a hyperspectral image . . . . . . 6 1.5 AVIRIS image of Cuprite, Nevada . . . . . . . . . . . . . . . . . . . . 7 1.6 Sliding window used in some algorithms . . . . . . . . . . . . . . . . 10 1.7 Performance comparison between xed and oating point . . . . . . . 11 1.8 Gantt diagramm of this work . . . . . . . . . . . . . . . . . . . . . . 12 2.1 Precision of calculations with dierent shift values . . . . . . . . . . . 19 3.1 A schematic of the top block . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Schematic of the implemented inverter . . . . . . . . . . . . . . . . . 23 3.3 Arithmeticpipeline............................ 24 3.4 Data depency in the pipeline . . . . . . . . . . . . . . . . . . . . . . . 25 3.5 Schematic of mean subtract. . . . . . . . . . . . . . . . . . . . . . . . 26 3.6 First proposed method for the matrix multiplication . . . . . . . . . . 27 3.7 Second proposed method for the matrix multiplication . . . . . . . . 27 3.8 Schematic of the dual matrix multiplier . . . . . . . . . . . . . . . . . 28 3.9 Schematic of the coordinate sorter . . . . . . . . . . . . . . . . . . . . 29 4.1 HYDICE sensor information and used dataset . . . . . . . . . . . . . 31 4.2 AVIRIS sensor information and used dataset . . . . . . . . . . . . . . 31 4.3 Spectral Angle Mapper algorithm . . . . . . . . . . . . . . . . . . . . 33 4.4 Results for the HYDICE dataset . . . . . . . . . . . . . . . . . . . . . 34 4.5 Results for the AVIRIS dataset . . . . . . . . . . . . . . . . . . . . . 35 iii List of Tables 2.1 DSP usage for multipliers with diferent widths . . . . . . . . . . . . . 18 3.1 Latencies for some arithmetic units . . . . . . . . . . . . . . . . . . . 24 4.1 Basic specications of the Virtex 690T FPGA . . . . . . . . . . . . . 30 4.2 FPGA resource usage overview . . . . . . . . . . . . . . . . . . . . . 36 4.3 Implementation latency overview . . . . . . . . . . . . . . . . . . . . 36 4.4 Final resource and time results for the implementation . . . . . . . . 37 iv Chapter 1 Introduction 1.1 Motivation and objectives Space exploration serves many purposes, the most obvious being gathering information about our planet and its surroundings. For this purpose, sensors capable of gathering information are created, such as antennas or telescopes that are used both from the Earth and sent aboard spaceships. One of those are hyperspectral cameras, which take pictures in hundreds of dierent bands. Their data allows to nd objects, detect materials or identify processes. As technology advances, these sensors evolve requiring appropriate processing solutions to interpret the data or compress it and send it to ground. The objective of this work is the implementation of one of these algorithms in a way that the processing in the aircraft is preferable to the transmission of the raw data. For this purpose, dierent algorithms will be evaluated and one will be chosen, more specically, the Reed Xiaoli algorithm. A rst implementation of the oating point algorithm will be done and its transformation to integer arithmetic will be studied. This step is necessary because the major impediment of these algorithms to be implemented in hardware is the high number and complexity of its operations. With the arithmetic well dened, its implementation will be adapted and a comparison of accuracy and performance between the rst oating point version and the integer version will be made. 1 1.2 Related work In this section we will review the state of the art on the use of FPGAs (eld programmable gate arrays) in space applications in general and implementations related to the algorithm developed in this work. In [15] a study is carried out on the current situation of the use of FPGA on board aircraft for hyperspectral analysis and implementations of two algorithms, ISRA and N-FINDR, are presented. The results show numerous advantages of FPGA over other types of solutions such as GPUs, such as its smaller size and weight and its resistance to radiation. In addition, their reconguration capabilities even after the system is launched are emphasized. In [11] an implementation of a target generation algorithm is presented. It uses an inverter based on the Gauss Jordan elimination method, same as the one that will be used in this work, and the results are presented on the same platform. In [9] a study of the same algorithm that will be implemented in this work is carried out on an FPGA. The results are positive with a reduction in calculation time compared to software-based solutions. However, this study was performed only on a oating point implementation which limited it to the processing of previously dimensionally reduced hyperspectral images. In [21] the processing of data on board is proposed with the aim of reducing network usage and accelerating data processing. One of the steps of this processing is also the RX algorithm that will be implemented here. In general, the previous works present good results on these systems and expectations to an increase of use, motivated by more and more advanced image capture systems that take the hardware to their limits. 2 Other methods Apart from the RX algorithm, there are other methods for detecting anomalies, although a notable number of them are based on RX. Subspace methods Subspace methods are global and apply principal component analysis (PCA) or singular value decomposition (SVD) to the hypercube. The rst PCA/SVD bands are supposed to be the background and are removed dierently by each method. Not removing enough bands means that the anomalies can be lost in the background noise and removing too many would lead the anomalies to disappear. The determination of the correct number of bands to be removed is still under study [6]. Subspace RX In this method, the RX algorithm is applied to a limited number of PCA bands. The rst components are discarded. RX after orthogonal subspace projection [16] In this method, the rst PCA/SVD components dene the subspace of the background and the data is projected onto an orthogonal space before applying XR. RX after partialling out the clutter subspace [14] In this method, clutter eect on a pixel is discarded component wise taking each of its spectral components as a linear combination of its high variance principal components. The RX algorithm is applied to the results. Complimentary subspace detector [8] This algorithm is not based on RX. In CSD, the principal components with higher variance are used to dene the subspace of the background and the other PCs to dene the subspace of the target. The pixel is then projected onto the two subspaces and the result is the dierence between them. Local methods In the local methods the background is derived from the neighboring pixels or surroundings of the pixel under test (PUT). Two windows are dened, a guard and an exterior, and the neighbors are the pixels that are between these two. Sometimes, for example in local RX, a third one is used where the covariance of the background is calculated in a window larger than the average of the background (see Figure 1.6). 9 Quasi-local RX [12] Quasi-local RX is a compromise between global RX and local RX, where the global covariance matrix is decomposed using eigenvector/eigenvalues. The eigenvectors are maintained, but the eigenvalues are replaced by the maximum local variance, resulting in lower detector scores in areas with high variance. Covariance window Mean window Guard window PUT Figure 1.6 : Sliding triple window used in the local AD methods. Segmentation based methods In complex scenes it is dicult to assume that the background will be dened by a normal distribution, so methods have been developed to separate it into dierent classes. Class-Conditional RX [13] In this method the image is rst segmented and the covariance and mean matrix calculated for each of these classes. Each pixel will correspond to the class in which its RX value is lower. There are more methods within those based on segmentation, from those using stochastic functions such as Method Based on Multivariate Normal Mixture Models to methods based on Self-organizing maps. 10 1.2.4 Comparing oating and xed point In computer systems, there are two numerical representations for real numbers, xed and oating point. They have dierent arithmetic, which gives them dierent range and resolution with the same number of bits. Therefore, there are certain applications or platforms more akin to one of them. For example, the ability of oating point numbers to contain in the same number of bits both very large and very small numbers and to adjust their resolution accordingly is very attractive from the programmer's point of view but the simplicity of xed point operations allow their use in small microcontrollers or to save resources in FPGAs (see Figure 1.7). Figure 1.7 : A FIR lter, originaly implemented as a single-precision oating-point lter and converted to xed-point. The xed-point design shows both resource reduction and latency improvements Due to the high complexity of hyperspectral image analysis, current FPGAs have little capacity to implement some of these algorithms. Therefore, one of the objectives of this work is to perform two parallel implementations, one in oating point and another one in xed point and compare their results. 11 1.3 Project plan First, an algorithm and a hardware platform are chosen. Desirably, this algorithm should be known and commonly found in the literature. This algorithm will be implemented in software and this implementation tested with real images against existing software such as ENVI or Spectral Python. The dierent algorithm stages will then be optimized for a hardware architecture and the transformation of the algorithm from oating point to integer or xed point logic will be studied. This transformation requires design decisions that sacrice precision with the goal of saving hardware resources, so the underlying hardware will have to be considered. Subsequently, a hardware validation of the design will be realized. Finally, a study will be carried out on the accuracy of the results obtained. Tasks 1 2 3 4 5 6 7 8 9 10 1 Research algorithms 2 Software impl. 3 Model hardware 4 Hardware impl. 5 Test and validation 6 Report results Figure 1.8 : Gantt diagram with an overview of months and tasks 12 Chapter 2 Software Model The RX algorithm has been chosen for this work because it is the benchmark for this kind of algorithms [6], [19], [16] and many existing algorithms basing on RX in some manner. In order to become familiar with the algorithm and create a platform where tests can be easily performed, a software implementation has been made as a rst step. The rst step is to divide the algorithm into simpler operations:  Calculate the mean, deviation and K covariance matrix of the image  Calculate K−1 that is, the inverse of the covariance matrix  Calculate δRX for each pixel in the image  Sort the results 13 2.0.1 Mean, deviation and covariancce matrix One of the bottlenecks in this type of FPGA-based systems is the input and output of data [20]. Since the calculations of mean, variance and covariance matrix need the original matrix -a cubefor this, it has been decided to calculate these in a CPU and transmit the results to the FPGA. In addition, the operations are relatively simple for a CPU. Next, the pseudocode of the three is presented. Algorithm 1 Pseudocode for the mean, deviation and covariance matrix 1: . bands : number of bands in the image 2: . pixels : number of pixels in the image 3: . A : the image in the form of a 2D matrix with pixels ×bands 4: .T is used to denote the transpose of a matrix 5: function mean ( A ) 6: for i←0 to bands −1 do 7: sum ←0 8: for j←0 to pixels −1 do 9: sum ←sum +A[i][j] 10: mean[i]←sum/pixels return mean 11: function deviation ( A, mean ) return AT−mean . This operation is a matrix subtraction 12: function covariance ( deviation ) return deviationT∗deviation/(pixels −1) The covariance, a band2 matrix is then transmitted to the FPGA to calculate its inverse. 14 2.0.2 Inverse Choosing the algorithm Before starting the implementation, several algorithms to perform the inverse have been studied. QR Factorization [4]: QR factorization breaks down the A matrix into the product of two matrices A=QR , with Q being an orthogonal matrix and R a superior triangular matrix. With this triangular matrix it becomes easy to calculate the inverse. However, although QR factorization can be eciently performed on a powerful matrix multiplication module or several modules that can be executed simultaneously, such as in a GPU, it does not take advantage of the capabilities provided by our system such as arbitrary width arithmetic units. Gauss Jordan elimination method [11]: The Gauss Jordan method dictates that if we have a A matrix that can be transformed into the identity matrix through elementary operations, these same operations transform the identity matrix into A−1 . Since it is possible to execute these elementary operations in an entire row at once and the operations between rows are independent, this method is easily parallelizable. Therefore, this was the chosen method. Generally speaking, the execution of the algorithm takes place in such a way that: 1. An identity matrix is generated 2. The same operations are performed on both matrices until the A matrix is transformed into the identity matrix 3. The result is in the matrix generated in the rst step As seen in the following pseudocode, these elementary operations are performed in 3 steps: 15 Algorithm 2 Pseudocode of the Gauss Jordan method 1: A : a square matrix with the size n∗n 2: A−1 : an identity matrix with the size n∗n 3: function inverse ( A ) 4: for i←0 to n−1 do . Forward elimination to build the upper triangular matrix 5: . row i acts as the pivot 6: if A[i][i]=0 then . If the later divisor is 0 7: for j←i+ 1 to n−1 do 8: if A[i][j]6= 0 then 9: A[i]←A[j], A[j]←A[i] 10: for j←0 to n−1 do 11: A−1[j]←A−1[j]−A−1[i]∗(A[j][i]/A[i][i]) 12: A[j]←A[j]−A[i]∗(A[j][i]/A[i][i]) . These two lines run in parallel 13: . After a complete iteration of the outer loop, the pivot contains the desired form 14: 15: 16: for i←n−1 to 0 do . Backward elimination to build a diagonal matrix 17: for j←i−1 to 0 do 18: A−1[j]←A−1[j]−A−1[i]∗(A[j][i]/A[i][i]) 19: A[j]←A[j]−A[i]∗(A[j][i]/A[i][i]) . These two lines run in parallel 20: 21: 22: for i←0 to n−1 do . Last division to build identity matrix 23: A−1[i]←A−1[i]∗(1/A[i][i]) 24: A[i]←A[i]∗(1/A[i][i]) 25: . There is no need to update the values in the starting matrix return A−1 16 2.0.3 Matrix multiplication With the inverse of the convariance matrix calculated in the FPGA, the real RX algorithm can continue to be performed, where a hyperspectral pixel is multiplied with this inverse. This calculation will require the CPU to transmit the original image, the original image next to the average, or the deviation already calculated to the FPGA. The RX algorithm dictates that: rx(x)=(x−µ)TK−1 N×N(x−µ) K−1 N×N being the inverse matrix with a dimension of N×N , (x−µ) being the deviation of a pixel a dimension of N×1 and (x−µ)T it's transpose, with a dimension of 1×N . In this step, the inverse matrix gets eectively multiplied by a unique pixel across all bands. Through row reduction, a single value for this pixel is recovered which can then be mapped to a 2d image. Since the operands are heterogeneous, use can be made of the associative law on matrix products that dictates that: A∗(B∗C) = (A∗B)∗C to optimize the operation. In the implementation section, a study on what alternative to use will be carried out, as this depends heavily on the hardware architecture. 2.1 Adapting the arithmetic As shown before in Figure 1.7, oating point arithmetic operations are very inecient compared to xed point operations. In the next chapter a transformation strategy will be dened and the steps followed to make this transition will be shown. Taking advantage that the RX algorithm only needs relative and not exact values, ie, the highest value found in the results will be the most anomalous (see 1.2.3) regardless of whether it exceeds a certain threshold or not, the xed-point representation may be simplied by ignoring the fractional part and keeping only the integer. 17 2.1.1 DSP Blocks The DSP blocks are prefabricated circuits found next to the FPGA logic and implement arithmetic operations. These blocks allow to perform operations more e- ciently and operate at higher frequencies than equivalent logic in the fabric, besides not occupying space in it. However, as they are prefabricated blocks, they do not oer the same exibility as the rest of the FPGA logic. This implies that the size of the operands and the cost of operation resources will not always be proportionally related. Therefore, dierent widths have been tested and their resource usage has been noted. The results of these tests for multiplication are in the following table: Width operand A Width operand B DSP48 Slices used 25 18 1 35 25 2 52 24 3 42 35 4 64 25 4 52 42 6 Table 2.1 : Depending on the width of the operands, the operation will use a dierent amount of resources (data obtained for signed multiplication in Vivado 2019.2 on Xilinx Series 7 chips) For subtraction operations up to 48bits, only one DSP is used, so they are not expected to require more than one block. The divisions cannot be implemented by this type of blocks, so they will use standard logic and can be used with an arbitrary width. With this, data for all the operations that are going to be carried out in the FPGA is recorded. 18 One of the operands is the row , which is used at two dierent times within the algorithm (see Figure 3.4). Since a second read on the BRAM cannot be done since it is busy performing the read of the next operations, this data is stored at the time of reading in a FIFO and will be fetched when needed. A−1[j]←A−1[j]−A[i]∗A[j][i]/A[i][i] Figure 3.4 : As can be seen in the algorithm, row j is used at two dierent times. Additionaly, it can be observed that in the rst step in which the upper triangular matrix is constructed, the algorithm requires a check on the pivot row and a possible exchange of rows. This is necessary because this value is later the dividend, so a 0 would cause a failure in the calculation. BRAM memory readings have a latency cycle, so reading an inappropriate value two cycles in a row -in the case of reading an inappropriate and then an appropriate value, there would be the possibility of performing an inplace row swap with the pivot and the row just after itwould not only add latency to the calculation but also increase the complexity of the module. Therefore, the dividend checks are done on the writes, recorded in a renaming table that will be checked at the time of reading. This ensures that the reads will always be valid for the calculation. In the case of the very rst division, this dividend comes directly from the CPU and the upper module control is responsible for reordering this row if necessary. This renaming table is located in registers, so it is possible to access it without any latency and as it only contains indexes, it does not overload the FPGA resources. In addition, this table is local, so the results have to be reordered in the RAM itself before leaving the inverse module. Since these swaps only occur in the calculation of the upper triangle, the lower triangle can be used to reorder them. The reordering system is very simple, the data enters the pipeline according to the order that exists in the renaming table and is written in its natural order. This implies that the results of this reordering will be correct as long as both rows that have been rotated are at the same time in the processing pipeline, which in the case of xed point is approximately 90 in size. Experimental results show that it is rarely necessary to rotate rows -although enough to recommend the inclusion of a method to deal with it-, and that these rotations rarely exceed one or two positions in the pipeline. In the transformation to integer arithmetic it was discovered that the calculation of the inverse is the one introducing more error in the nal results of the algorithm, 25 therefore an exhaustive study on how to minimize it has been made. For this it has been necessary to reduce the values in which the limited precision produced overows and increase small values to give them more weight in the operations. By performing the operations of identity generation, upper triangle, lower triangle and diagonal independently, it has been possible to place dierent shift values and further rene the resolution of the algorithm. These operations are performed by the shift process. It should also be said that there is an error in the generation of Xilinx dividers. When you enter numbers near the precision limit you lose control of the sign. Therefore, div_x , a process that converts all the entered operands into positives and saves their position in a pipeline has been placed before the division. When the results are produced, the tag is checked in the pipeline, the negative is calculated and replaced if necessary. 3.2.3 Mean subtract Figure 3.5 : Schematic of mean subtract. The mean subtract module receives the calculated average and the original pixels of the image and subtracts them. This calculation is the deviation and although it had already been calculated by the CPU, it is possible that the latter discards the data to free up space. The calculation of the average is required because it is assumed that its size being much smaller, the CPU can keep it in memory. In case the deviation can be received directly from the CPU, this module can be simply deleted. The module receives the elements from the upper module that reads them from the same FIFO. The rst elements are the average and are stored in a BRAM that is treated as a circular buer, and the following elements are directly subtracted and returned to the upper module again. 26 3.2.4 Matrix multiplication As noted above (see 2.0.3), this calculation can be implemented in two dierent ways. Following will be a comparison of the rst required multiplication in both methods: (x−µ)TK−1 N×N : A row from the inverse and the whole column of the deviation get read, each element multiplied with its correspondent and all products added together. If stalls were to be avoided, this sum would need to be computed every cycle, which can easily be achieved with an adder tree.   a b c d e f g h i  ∗  1 2 3  =  1∗a+ 2 ∗b+ 3 ∗c 1∗d+ 2 ∗e+ 3 ∗f 1∗g+ 2 ∗h+ 3 ∗i   Figure 3.6 : First proposed method for the computation of a single pixel: red, blue and green represent data processed in the rst, second and third cycles respectively. Note that the entire deviation data of that pixel gets used every cycle K−1 N×N(x−µ) : The inverse gets also read row by row, but the deviation matrix only by elements. Each element of the rst row of the matrix gets multiplied with the rst element of the deviation, the result accumulated, and continued with the next pair row/element. This goes for N cycles, that is, a whole inverse matrix and a whole pixel in the deviation matrix. The result is N accumulated values which get ushed every N cycles, which ends up being the same throughput as the former method.  123 ∗  a b c d e f g h i  = 1*a + 2*d + 3*g 1*b + 2*e + 3*h 1*c + 2*f + 3*i  Figure 3.7 : Second proposed method: red, blue and green represent data processed in the rst, second and third cycles respectively. Here only an element of the deviation data gets accessed each cycle. While both methods have equivalent cost in time -the former has the added latency of the adder tree, the latter the latency of the accumulatorsand also similar 27 cost in DSP usage, data input by row is less taxing on the CPU and its FIFO structure can be reused for the second multiplication. Henceforth, the second approach was chosen. Figure 3.8 : Schematic of the dual matrix multiplier The second multiplication is similar in both steps, a 1×N by N×1 multiplication. One operand comes every cycle and each N cycles all products get added together. This sum is realized through an accumulator. The module contains three subprocesses:  rst_mac reads the inverse and performs its multiplication with the received deviation. This deviation is also stored in a FIFO. The products are then accumulated till a whole pixel has been computed.  second_mac stores the results of rst_mac in registers and performs the multiplication with data from the FIFO, with the result being accumulated. Every cycle, the registers are shifted so a new multiplication is done.  write_proc controls the writing of the results from second_mac to the sorter and computes the coords. 28 3.2.5 Coordinate sorter Figure 3.9 : Schematic of the coordinate sorter This module receives a value and a pair of coordinates every bands cycles. These values are written in a BRAM memory that acts as an ordered list. Each entered value is compared with the head of the list, the highest value is saved and the other is saved in a temporary variable, compared with the second value in the list, and so on. Since a value is received each bands , the maximum number of possible values to be stored in this list is also bands . The rest of the values are discarded. When the last value has been introduced, the module communicates the highest pixels, that is, the most anomalous ones, to the superior module so that they are communicated to the CPU. 29 Chapter 4 Results 4.1 Recongurable platform The architecture described in the previous section has been implemented using the VHDL language. To test its correct operation, the Vivado environment and the Virtex 690T FPGA 4.1 have been used. Although it is not a radiation-protected FPGA, the architecture of this algorithm is scalable so that its adaptation to other sensor sizes or FPGAs should not be a problem. 4.2 Hyperpectral image datasets Two hyperspectral images have been used for the work, one taken by the HYDICE sensor and the other by the AVIRIS sensor. Both images are commonly used as reference in hyperspectral applications. In the eld of the scene captured by HYDICE, 15 panels of dierent sizes were placed on a eld in a 3 x 5 meter conguration. The following images show a false color image with the bands 50 , 37 , and 17 as red, green and blue respectively [10] and the location of the panels as detected by software. Part number Slices Logic cells Flip-Flops BRAM DSP Slices XCE7VX690T 108,300 693,120 866,400 1,470Kb 3,600 Table 4.1 : Basic specications of the Virtex 690T FPGA 30 Horizontal res 64 Vertical res 64 Bands 210 (169*) Spectral res 400 - 2500nm Spatial res 1,56meter/pixel Size 1,6 MB Figure 4.1 : HYDICE sensor information [3], image in false color and anomaly map. *In this work, images and results are reported for 169 bands The image taken by the AVIRIS sensor was taken on September 16, 2001, ve days after the terrorist attacks that brought down the WTC towers and its surrounding buildings. The spatial resolution of this image is very high because a very low altitude ight was performed. Along with the false color image, an image with the anomalies detected by software is provided. To ease the recognition of these anomalies, a circle has been painted around each group of anomalues for the rst 30. Horizontal res 614 Vertical res 512 Bands 224 Spectral res 360 - 2500nm Spatial res 1,7meter/pixel Size 140 MB Figure 4.2 : AVIRIS sensor information [2], image in false color and anomaly map 31 4.3 Adequacy of approximation 4.3.1 Floating point The results provided by the oating point version of this system are the same as those provided by the equivalent software version, so it can be considered valid. 4.3.2 Fixed point The xed point version of the system is an approach to the oating point system with the intention of maintaining the highest possible accuracy with limited use of resources. Therefore, the results are dierent and an assessment of their accuracy must be made. For this purpose, three metrics have been used. In the rst and simplest, it has been veried that the number of detected anomalies can be found in the rst x positions in both results, regardless of the order. Extending the previous strategy, it has been checked for the non-coinciding elements, if any neighbor has been detected in the environment. A neighbor is dened as an adjacent pixel, both in a straight line and in diagonal. Finally, the rst strategy has been extended again, this time it has been checked if for the non-coincidences an anomaly with a spectral similarity of less than 5 degrees has been found. Spectral similarity The Spectral Angle Mapper (SAM) is a physics-based spectral classication that uses an n-D angle to match pixels with reference spectra. The algorithm (see section 4.3.2) determines the spectral similarity between two spectra by calculating the angle between the spectra and treating them as vectors in a space with a dimensionality equal to the number of bands. This technique, when used in calibrated reectance data, is relatively insensitive to the eects of illumination and albedo. 32 α= cos−1     nb P i=1 tiri nb P i=1 t2 i1/2nb P i=1 r2 i1/2      α = spectral angle between vectors nb = number of spectral bands t = target pixel r = reference pixel Figure 4.3 : Spectral Angle Mapper algorithm 4.3.3 Results for HYDICE Every pixel is found in the same order in reference and simulated calculations. Also, since those pixels are found equally, it is obvious that they are also found as neighboring pixels and are spectrally similar. Since it is an image taken to test the sensor, the targets are big and clear, and the image is quite easy to analyze. 4.3.4 Results for AVIRIS The WTC image is more complex and the detection of exact pixels and neighbors quickly falls. However, the spectral signatures match for the most part, especially considering the number of real hot spots found in the image. Therefore, the results can also be taken as successful. 33 Figure 4.4 : A diagram showing similarities between referenced and achieved results and reference and achieved results mapped as in the original image for the HYDICE dataset 34 [19] Dalton Rosario. A semiparametric model for hyperspectral anomaly detection, January 2012. [20] Qingshan Tang, Habib Mehrez, and Matthieu Tuna. Multi-FPGA prototyping board issue: the FPGA I/O bottleneck. In 2014 International Conference on Embedded Computer Systems: Architectures, Modeling, and Simulation (SAMOS XIV) , pages 207214, July 2014. [21] James Theiler, Bernard Foy, Claira Sa, and Steven P. Love. Onboard CubeSat data processing for hyperspectral detection of chemical plumes. In David W. Messinger and Miguel Velez-Reyes, editors, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XXIV , page 5, Orlando, United States, May 2018. SPIE. [22] Xin Zhou, Norihiro Tomagou, Yasuaki Ito, and Koji Nakano. Ecient Hough Transform on the FPGA using DSP Slices and Block RAMs. In 2013 IEEE International Symposium on Parallel Distributed Processing, Workshops and Phd Forum , pages 771778, May 2013. 41