Hyperspectral image compression onboard next-generation satellites: implementation solutions on GPU and FPGAs

HYPERSPECTRAL IMAGE COMPRESSION ONBOARD NEXT-GENERATION SATELLITES: IMPLEMENTATION SOLUTIONS ON GPU AND FPGAS María Lucana Santos Falcón TESIS DOCTORAL Tesis Doctoral María Lucana Santos Falcón Las Palmas de Gran Canaria, Julio 2014 HYPERSPECTRAL IMAGE COMPRESSION ONBOARD NEXT-GENERATION SATELLITES: IMPLEMENTATION SOLUTIONS ON GPUs AND FPGAs UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA Instituto Universitario de Microelectrónica Aplicada DIVISIÓN DE DISEÑO DE SISTEMAS INTEGRADOS TESIS DOCTORAL Hyperspectral image compression onboard next-generation satellites: implementation solutions on GPUs and FPGAs María Lucana Santos Falcón Acknowledgements I would not have been able to complete this journey without the aid and support of countless people over the past four years. Foremost, I would like to express my gratitude to my supervisors, Prof. Jos´e L´opez and Prof. Roberto Sarmiento, who have been greatly supportive and have guided me during this research and while writing this Thesis, offering constructive comments and warm encouragement. Over the years, I have received funding from several entities, which have supported me while I completed my PhD. I would like to thank Thales Alenia Space Espa˜na S.A., the European Network of Excellence on High Performance and Embedded Architecture and Compilation (HIPEAC), the Institute for Applied Microelectronics (IUMA) and Ayuntamiento de Las Palmas de Gran Canaria for their financial support. I highly appreciate the feedback offered by Luis Berrojo and Javier Moreno, from Thales Alenia Space Espa˜na S.A., who have made an important contribution to the technical quality of this Thesis. I also thank the European Space Agency for offering me the opportunity to make a fruitful research stay. In particular I would like to thank Raffaele Vitulli and Luca Fossati for their generous support and for sharing their immense knowledge. I am also grateful to the Group on Interactive Coding of Images (GICI) from Universitat Aut`onoma de Barcelona and to Enrico Magli, from Politecnico di Torino for kindly offering their help. I’d like to thank also my fellow labmates, for all the stimulating discussions, the fun we have had, the coffees and their patience. I could not have imagined a better work environment for pursuing the PhD. Last, but not least, I am deeply grateful to my parents for generously offering me support and the education that has made it possible for me to get here. Thanks also to my sister, Isabel, the rest of my family and to Lola and Ricardo, for instilling me confidence and believing in me. vii Contents Abstract i Resumen iii Acknowledgements vii List of Figures xv List of Tables xix Abbreviations xxiii Symbols xxviii 1 Introduction 1 1.1 Outline .............................. 2 1.2 Preliminary concepts . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Multispectral and hyperspectral images . . . . . . . . 4 1.2.2 Instruments and sensors for hyperspectral data collection 6 1.2.3 Applications of hyperspectral images . . . . . . . . . . 8 1.2.4 Hyperspectral image compression . . . . . . . . . . . . 11 1.2.4.1 Lossless versus lossy compression . . . . . . . 14 1.3 Motivation of research . . . . . . . . . . . . . . . . . . . . . . 15 1.3.1 Importance of on-board hyperspectral image compression ............................ 16 ix x Contents 1.3.2 Limitations and difficulties of the on-board hardware . 18 1.4 Researchgoals .......................... 21 1.5 Organization of this document . . . . . . . . . . . . . . . . . 22 2 On-board hyperspectral image compression algorithms and hardware implementations 25 2.1 Outline .............................. 26 2.2 Algorithms for on-board hyperspectral image compression . . 27 2.2.1 Requirements and limitations of an on-board hyperspectral image compression algorithm . . . . . . . . . 30 2.2.2 Transform-based compression algorithms for hyperspectralimages ........................ 31 2.2.3 Prediction-based compression algorithms for hyperspectralimages ........................ 34 2.2.4 Recent research on hyperspectral image compression algorithms......................... 38 2.2.5 CCSDS Standard algorithms for satellite data compression.......................... 39 2.3 Physical implementations for on-board compression of hyperspectralimages.......................... 41 2.3.1 On-board hardware technology requirements . . . . . 43 2.3.2 Software implementations . . . . . . . . . . . . . . . . 47 2.3.2.1 Implementations on general-purpose central processing units (CPUs) . . . . . . . . . . . . 47 2.3.2.2 Implementations on digital signal processors (DSPs) ..................... 49 2.3.2.3 Implementations on graphics processing units (GPUs)..................... 51 2.3.3 Hardware implementations . . . . . . . . . . . . . . . 53 2.3.3.1 Hardware design flow . . . . . . . . . . . . . 53 2.3.3.2 Implementations on application-specific integrated circuits (ASICs) . . . . . . . . . . . . 54 2.3.3.3 Implementations on field-programmable gate arrays (FPGAs) . . . . . . . . . . . . . . . . 56 3 Implementation of a lossy compression algorithm for hyperspectral images on a GPU 65 3.1 Outline .............................. 66 Contents xi 3.2 LCE algorithm description . . . . . . . . . . . . . . . . . . . . 67 3.2.1 Prediction......................... 69 3.2.2 Rate-distortion optimization . . . . . . . . . . . . . . 71 3.2.3 Quantization and mapping . . . . . . . . . . . . . . . 72 3.2.4 Entropy coding . . . . . . . . . . . . . . . . . . . . . . 72 3.2.5 Fileformat ........................ 73 3.2.6 LCE compression efficiency . . . . . . . . . . . . . . . 74 3.3 Software implementation of the Lossy Compression for Exomars (LCE) algorithm . . . . . . . . . . . . . . . . . . . . . . 74 3.3.1 Generation of the compressed file . . . . . . . . . . . . 76 3.3.2 Configuration parameters . . . . . . . . . . . . . . . . 77 3.4 GPU architecture and NVidia CUDA . . . . . . . . . . . . . 78 3.4.1 CUDA abstractions . . . . . . . . . . . . . . . . . . . 79 3.4.2 CUDA memory spaces . . . . . . . . . . . . . . . . . . 81 3.4.3 Nvidia TESLA C2075 GPU . . . . . . . . . . . . . . . 81 3.5 Parallelization of the LCE compressor with Computer Unified Device Architecture (CUDA) . . . . . . . . . . . . . . . . . . 85 3.5.1 Allocation of the image data in the GPU . . . . . . . 86 3.5.2 Prediction, quantization and mapping . . . . . . . . . 86 3.5.3 Entropy coding . . . . . . . . . . . . . . . . . . . . . . 88 3.5.4 Bitpacking ........................ 94 3.6 Parallelization of the LCE decompressor . . . . . . . . . . . . 101 3.6.1 Preliminary considerations . . . . . . . . . . . . . . . 102 3.6.2 Header design . . . . . . . . . . . . . . . . . . . . . . . 103 3.6.3 Decoding the blocks . . . . . . . . . . . . . . . . . . . 107 3.6.4 Inverse quantization and prediction . . . . . . . . . . . 108 3.7 Experimental results . . . . . . . . . . . . . . . . . . . . . . . 108 3.7.1 Validation......................... 109 3.7.2 Impact of adding a header . . . . . . . . . . . . . . . . 110 3.7.3 Profiling.......................... 111 3.7.4 Speedup.......................... 113 3.7.5 Throughput........................ 115 3.7.6 Effect of the configuration parameters in the performance of the GPU implementation of the LCE compressor........................... 118 xii Contents 4 Implementation of a lossy compression algorithm for hyperspectral images on an FPGA 121 4.1 Outline .............................. 122 4.2 CatapultC design flow . . . . . . . . . . . . . . . . . . . . . . 124 4.3 Adapting the C language source code of the LCE algorithm forCatapultC........................... 128 4.3.1 Identification of the top function and inputs and outputs of the design . . . . . . . . . . . . . . . . . . . . 129 4.3.2 Configuration parameters . . . . . . . . . . . . . . . . 132 4.3.3 Reducing the complexity of the mathematical operations to calculate the gain factor α........... 133 4.3.4 Loop optimization . . . . . . . . . . . . . . . . . . . . 137 4.4 Results of the FPGA implementation of the LCE algorithm withCatapultC.......................... 141 4.4.1 Manual scheduling of the design . . . . . . . . . . . . 145 4.4.2 Implementation of the LCE algorithm using a modular approach ......................... 147 4.4.3 Comparison with the FPGA implementation of a nearlossless algorithm . . . . . . . . . . . . . . . . . . . . . 148 4.5 Performance comparison: FPGA, GPU, CPU . . . . . . . . . 149 5 Implementation of the CCSDS standard for lossless hyperspectral image compression on a space-qualified FPGA 155 5.1 Outline .............................. 156 5.2 The CCSDS 123 standard for lossless multispectral and hyperspectral image compression overview . . . . . . . . . . . . 157 5.2.1 Prediction......................... 159 5.2.2 Entropy coding . . . . . . . . . . . . . . . . . . . . . . 163 5.3 Design methodology . . . . . . . . . . . . . . . . . . . . . . . 163 5.4 Impact of the user-defined parameters in the compression efficiency .............................. 165 5.5 Architectural design considerations . . . . . . . . . . . . . . . 169 5.5.1 Encoding order . . . . . . . . . . . . . . . . . . . . . . 171 5.5.2 Local sum mode and prediction mode . . . . . . . . . 172 5.5.3 Number of bands for prediction . . . . . . . . . . . . . 173 5.5.4 Hardware complexity estimation . . . . . . . . . . . . 174 5.5.4.1 Comparison and complexity estimation . . . 177 5.6 HyLoC Hardware architecture description . . . . . . . . . . . 181 Contents xiii 5.6.1 HyLoC verification and validation . . . . . . . . . . . 185 5.7 Experimental results . . . . . . . . . . . . . . . . . . . . . . . 189 5.8 Comparison of hardware technologies for the implementation of hyperspectral image compression algorithms . . . . . . . . 191 5.8.1 Implementations on GPUs . . . . . . . . . . . . . . . . 192 5.8.2 Implementations on FPGA . . . . . . . . . . . . . . . 193 5.8.2.1 Comparison with state-of-the-art FPGA implementations of hyperspectral compression algorithms ................... 196 5.8.3 Implementations on space-qualified FPGAs . . . . . . 198 6 Conclusions 201 6.1 Further research work . . . . . . . . . . . . . . . . . . . . . . 207 A Sinopsis en espa˜nol 211 A.1 Introducci´on............................ 212 A.2 Objetivos y metodolg´ıa de trabajo . . . . . . . . . . . . . . . 215 A.3 Compresi´on con p´erdidas en GPU y FPGA . . . . . . . . . . 217 A.3.1 Implementaci´on del algoritmo LCE en una GPU . . . 219 A.3.2 Implementaci´on del algoritmo LCE en una FPGA . . 220 A.4 Compresi´on sin p´erdidas en FPGA . . . . . . . . . . . . . . . 224 A.4.1 Algoritmo CCSDS 123 . . . . . . . . . . . . . . . . . . 225 A.4.2 Implementaci´on del CCSDS 123 sobre una FPGA cualificada para el espacio . . . . . . . . . . . . . . . . . . 226 A.5 Conclusiones ........................... 227 B Publications 229 B.1 Journals.............................. 230 B.2 International Conferences . . . . . . . . . . . . . . . . . . . . 230 References 233 List of Figures 1.1 Hyperspectral data cube. . . . . . . . . . . . . . . . . . . . . 4 1.2 Electromagnetic spectrum. . . . . . . . . . . . . . . . . . . . . 5 1.3 Remote sensing scanners. a) Whiskbroom b) Pushbroom . . 7 1.4 (a) Band-sequential order (b) band-interleaved order . . . . . 9 1.5 Example of the generation of the codeword with Golomb codes 13 2.1 Prediction-based compression scheme. . . . . . . . . . . . . . 34 2.2 CCSDS 123 Recommendation for lossless multiand hyperspectral image compression. . . . . . . . . . . . . . . . . . . . 41 2.3 Radiation effects on a MOSFET transistor.(a) Normal operation. (b) Post irradiation. . . . . . . . . . . . . . . . . . . . 45 2.4 LEON3 spacecraft controller on a chip. . . . . . . . . . . . . . 49 2.5 Simplified CPU and GPU architecture comparison. . . . . . . 52 2.6 CWICOM compression ASIC. . . . . . . . . . . . . . . . . . . 55 2.7 Basic elements of an FPGA. . . . . . . . . . . . . . . . . . . . 56 3.1 Division of the hyperspectral cube into blocks and notation. . 69 3.2 Prediction neighbourhood of the LCE algorithm: a) first band b)allotherbands. ........................ 70 3.3 File format of the LCE compressed data. . . . . . . . . . . . . 74 3.4 Rate-distortion curves for AVIRIS when compressed with LCE and other algorithms of the state-of-the-art [33] . . . . . . . . 75 3.5 Flowchart of the LCE algorithm . . . . . . . . . . . . . . . . 76 3.6 Pseudo-code of the main function of the LCE algorithm implementation in C language. . . . . . . . . . . . . . . . . . . . 77 3.7 CUDA abstractions: threads, blocks and grid. . . . . . . . . . 80 3.8 CUDA memory spaces . . . . . . . . . . . . . . . . . . . . . . 81 xv Abbreviations 2D two-dimensional 3D three-dimensional ASIC application-specific integrated circuit BIP band interleaved by pixel BPE bit plane encoder CALIC Context-based, Adaptive, Lossless Image Codec CCD charge-coupled device CCSDS Consultative Committee for Space Data Systems CMOS complementary metal-oxide semiconductor CNES Centre National d’Etudes Spatiales COTS commercial off-the-shelf CPU central processing unit CUDA Computer Unified Device Architecture xxiii xxiv Abbreviations DCT Discrete Cosine Transform DPCM differential pulse code modulation DSP digital signal processor DWT Discrete Wavelet Transform EDAC error detection and correction EEPROM electrically erasable programmable read only memory EPROM erasable programmable read only memory ESA European Space Agency FAPEC Fully Adaptive Prediction Error Coder FET field effect transistor FL Fast-Lossless FLOPS floating-point operations per second FP7 The Seventh Framework Programme FPGA field-programmable gate array FSM finite-state machine GDDR Graphics Double Data Rate GPU graphics processing unit HLS high-level synthesis ITI Innovation Triangle Initiative Abbreviations xxv JPL Jet Propulsion Laboratory KLT Kahrunen-Lo`eve Transform LCE Lossy Compression for Exomars LUT lookup-table LWIR long wavelength infrared MAE maximum absolute error MSE mean-squared error MSRE mean square-root error MWIR medium wavelength infrared NASA National Aeronautics and Space Administration NIR near infrared NOAA National Oceanic and Atmospheric Administration OpenMP Open Multi-Processing PC principal component PCA Principal Component Analysis PLB Processor Local Bus PSNR peak signal-to-noise ratio RAM random access memory RD rate-distortion xxvi Abbreviations RD rate-distortion RTL register-transfer level SDRAM synchronous dynamic random access memory SEE single event effects SEU single event upset SLSQ Spectrum-oriented Least Squares SM streaming multiprocessor SMID single instruction multiple data SoC System-on-Chip SPECK set partitioned embedded block SPIHT set partitioning in hierarchical trees SRAM static random access memory SWIR short wavelength infrared TD Tucker Decomposition TDP thermal design power TID total ionizing dose TIR thermal infrared TMR triple modular redundancy UTQ uniform-threshold quantizer Abbreviations xxvii VNIR visible and near infrared VSWIR visible shortwave infrared Symbols xrow index yline index zband index nrow index within a block mline index within a block simage sample ˆspredicted sample ˜sreconstructed sample αleast-square estimator µaverage value α0quantized value of α µ0quantized value of µ ez,y,x prediction error kGolomb-code parameter Nx number of rows Ny number of lines Nz number of bands Nnumber of pixels in a squared spatial block xxix xxx Symbols hb horizontal block index vb vertical block index Bz,hb,vb N×Nblock with spatial coordinates (hb, vb) in band z Nhb number of horizontal blocks Nvb number of vertical blocks NB total number of blocks iblock index i=hb +vb ×Nhb jsample index within a block of N×Nsamples j=m+n×N qjnumber of bits taken by the j-th codeword of a block Qjbit position of the j-th codeword in the compressed block libits left unused in the last 32-bits word of the i-th compressed block Liprefix-sum of li piposition of the last codeword within block i Piword position where compressed block istarts in the compressed stream CH chunk NC number of chunks Kfirst element of the last chunk dxeceil xto the nearest integer bxcfloor xto the nearest integer sh leftinumber of bits that block ihas to be shifted to the left sh rightinumber of bits that block ihas to be shifted to the right Dedicado a Isabel, Lucana, Juan y Ricardo xxxi 6 Chapter 1. Introduction Although multispectral and hyperspectral images have different features, for the sake of simplicity, in the rest of this document we will use the term “hyperspectral” to refer to both. 1.2.2 Instruments and sensors for hyperspectral data collection The data collection in hyperspectral systems is a four dimensional problem, consisting of two spatial, one spectral and one time dimension, i.e. the collection of a hyperspectral cube requires scanning, in a time interval, in the spectral or spatial domain. In a remote sensing instrument, the light reflected by objects is captured by photodectector arrays, in which each element records one pixel of the image. Remote sensors, as opposed to photographic, acquire data using scanning systems that sweep over the terrain to build up and produce a two-dimensional image of the surface. The hyperspectral scanners are of two main types: whiskbroom and pushbroom, as depicted in Figure 1.3. Whiskbroom scanners or across track scanners reflect light into a single dectector using a mirror which moves back and forth to collect measurements from one pixel in the image at a time. All LANDSAT sensors prior to LANDSAT 8 use the whiskbroom design [1]. Other examples of whiskbroom imagers are the ones from the National Oceanic and Atmospheric Administration (NOAA): the Advanced Very High Resolution Radiometer (AVHRR) and the Geostationary Operational Environmental Satellite (GOES). On the other hand, pushbroom scanners, or along track scanners, use a line of detectors arranged perpendicular to the flight direction of the spacecraft. The image is collected one line at the time as the spacecraft flies forward. Instruments that use the pushbroom design include the ones on-board LANDSAT 8 [1], SPOT [2], IRS [3], PLEIADES [8] and PROBA [9]. The selection of the type of scanner depends mainly on the purpose and specific requirements of the Chapter 1. Introduction 7 remote sensing mission. Whiskbroom scanners are mechanically more complex but simpler from the optical point of view, while pushbroom imagers consist of a simpler mechanical system but more complex optics. The swath width is longer for whiskbroom imagers, when compared with pushbroom. Moreover pushbroom sensors must have the detectors perfectly calibrated to avoid stripes in the data caused by the varying sensitivity in the along track direction. Rotating mirror n bands Prism Swath width Scan direction Flight direction Swath width Flight direction Optics Array of detectors Dispersing element n bands (a) (b) Figure 1.3: Remote sensing scanners. a) Whiskbroom b) Pushbroom As it was already mentioned, hyperspectral imagers have to collect information of each pixel at different wavelengths. There are several types of devices for the spectral division or selection, which can be divided in three main classes: dispersive spectrometers; Fourier Transform interferometers; and narrow band tunable filters. Dispersive spectrometers use grating or a prism coupled with a two-dimensional array of detectors, in such a way that a spectral image is produced with the spatial information along one axis and 8 Chapter 1. Introduction the spectral information along the other. Fourier Transform interferometers split the radiation into two beams, introducing a controlled phase shift, and recombining them. The wavefronts of the beams on recombination interfere by the principle of superposition and the combined beam is focused on a detector. Finally, the narrow band tunable filters pass radiation through a very narrow bandpass or spectral bin, which can be spectrally tuned over different wavelengths, usually in a very short time. Although most past and current hyperspectral sensors have been airborne, many new space-based hyperspectral sensors have been proposed recently. These sensors have become increasingly important, because they are able to achieve near global coverage repeated at regular intervals of time, providing more imagery than their airborne counterpart. Once the image samples are detected they are converted into digital values, which can be arranged in different forms, commonly known as bandsequential or band-interleaved formats. In band sequential order, the samples are stored in raster order, band by band, i.e. all the samples in a specific band are stored before the storing of the samples in the next band start. On the other hand, in band interleaved format, the information is stored in such a way that all the spectral information of a specific pixel or line of pixels is stored before continuing with the next pixel or line of pixels. This is better illustrated in Figure 1.4. 1.2.3 Applications of hyperspectral images Hyperspectral technologies were developed for science and research purposes, but have progressed to allow the development of a diversity of commercial applications. The captured three-dimensional data cubes are processed after they are collected, so that the most relevant information can be extracted from them [10]. This information is extremely useful in many diverse fields, Chapter 1. Introduction 9 Ny Nx Nz bandn band0 band1 band2 band3 (b) Band-Interleaved order (a) Band-sequential order Figure 1.4: (a) Band-sequential order (b) band-interleaved order including mineralogy, Earth monitoring, surveillance and medicine. Some of the applications of hyperspectral data analysis are listed next. Atmospheric characterization and climate research Geologic mapping Environmental monitoring Vegetation analysis, food safety Monitoring of coastal environment Urban grow analysis Biological and chemical detection Non invasive diagnosis in cancer detection Surveillance 10 Chapter 1. Introduction Detection of weapons of mass destruction Detection of landmines Some of the main techniques utilized to extract useful information from hyperspectral data are summarized next. Dimensionality reduction. The dimensionality of the input hyperspectral image is reduced, in order to facilitate the subsequent processing of the scene [11]. Hyperspectral unmixing. The signal captured by a hyperspectral sensor at a given band and from a given pixel is a mixture of the reflectances of the different materials located in the respective pixel area. When mixing occurs, it is not any more possible to determine which materials are present in the pixels directly from the measured spectral vectors. The ability to discriminate materials can be recovered by applying hyperspectral unmixing techniques. An extensive amount of research work has been devoted to hyperspectral unmixing [12, 13]. Classification. Hyperspectral image classification has been a very active area of research in recent years [14]. Given a set of observations, the goal of classification is to assign a unique label to each pixel vector, so that it is well-defined by a given class. Hyperspectral target detection. Hyperspectral imagery has been used in reconnaissance and surveillance applications where targets of interest are detected and identified [15]. In the process of detecting a target, first the anomalies in the image are extracted [16, 17]. Then, the targets can be identified by their spectral signature, by comparing it with the data available in a spectral library [18] or from a set of training data. Chapter 1. Introduction 11 1.2.4 Hyperspectral image compression Hyperspectral images acquired by aircrafts or satellites represent a high amount of data. For instance, a single AVIRIS image occupies 134 Mbytes of data storage. In the specific case of remote sensors placed on satellites, the computational power, storage and downlink bandwidth are limited, therefore applying compression techniques have become an effective and cheap solution in order to cope with these restrictions and still allow the hyperspectral imagers to produce images at a high data rate. Compression of hyperspectral images is effective because the pixel values of neighbouring locations and wavelengths are highly correlated. Removing the correlation allows for reducing the data volume. Hyperspectral images have similar characteristics to natural photographic images or video, and consequently their size can be reduced with compression tools which were developed for image or video [19, 20]. The fundamentals of data compression are based on representing the necessary information with the smallest possible amount of bits. In general, compression can be lossless or lossy. Lossless techniques make it possible to recover all the original information after decompression, while lossy methods permit to recover the original data with some losses of information, enabling higher compression ratios. Compression techniques take advantage of several facts. The spatial or spectral redundancy in the images makes it possible to deduct a pixel value by using information of neighbouring pixels. Decorrelation methods like prediction o transformation can be applied, in such a way that it is no longer necessary to transmit the whole pixel information, but just the necessary information for the decoder to be able to calculate - or approximate if the technique is lossy - the value of the original pixel. In addition, compression methods exploit the statistical redundancy in the data, exploring the 12 Chapter 1. Introduction probability of the symbols in such a way that long codewords are used to represent symbols with low probability and short codewords are utilized to represent the most frequent symbols. These codes try to reduce the redundancy present in a source of information, and represent it with fewer bits which carry more information, minimizing the average length of the messages according to a particular assumed probability model, which is known as entropy encoding. Among the most common entropy encoding techniques we can find: Huffman codes. In Huffman codes [21], the source is encoded using a table of codes of variable length, which has been derived based on the estimated probability of occurrence for each possible value of the source symbol. Arithmetic codes. These type of codes convert strings of data into single floating point numbers between 0 and 1. They establish a model of the entire data set and find the occurrences of sequences of symbols that can be expressed in the form of a single number with high precision. Universal codes. Universal codes are prefix codes which map positive integers into binary codewords, ensuring that the length of the resulting codeword is within a constant factor of the expected lengths that an optimal code would have assigned. Golomb codes. These codes are prefix codes, used when the distribution of the source data is geometric. Golomb codes [22] first find the quotient and the remainder of the division of the source and a tunable parameter. The codeword consists of the quotient of the division expressed in unary notation, followed by a stop bit and the remainder in truncated binary notation, as showed in Figure 1.5. These codes are known as Rice codes [23] when the tunable parameter (the divisor) is a power of two. Golomb and Rice codes are widely used for lossless Chapter 1. Introduction 13 image and video compression, as well as for satellite data compression. In fact, it is employed in the Consultative Committee for Space Data Systems (CCSDS) standard for universal lossless satellite data compression, CCSDS 121 [24]; and in the standard for multispectral and hyperspectral image compression CCSDS 123 [25]. N Input sample q = floor (N/m) r = N-q*M Codeword q = floor (N/m) q ones r least significant bits of N Figure 1.5: Example of the generation of the codeword with Golomb codes Regardless the selected type of entropy encoding, there is a limit which determines the smallest possible expected number of bits needed to encode an event, known as the Shannon limit [26]. Shannon limit establishes that, given a set of mutually distinct events e1, e2, e3...en, and the probability distribution Pof the events, the smallest possible expected number of bits needed to encode an event is the entropy of P, denoted by: H(P) = n X k=1 −p{ek}log2p{ek}(1.1) In the former equation, p{ek}is the probability that event {ek}occurs. An optimal code outputs −log2pbits to encode an event whose probability of occurrence is p. It is important to assess how well a specific compression method performs when it is employed, since the compression efficiency will depend on the compression technique as well as on the distribution of the image data. The 14 Chapter 1. Introduction compression ratio is calculated and expressed commonly in two different ways: as the relationship between the size in bits of the raw image and the compressed image; or as the number of bits necessary to represent a pixel after compression, also known as bits per pixel. The formulas that can be employed to estimate how much a compression method reduces the data volume are shown next: Compression ratio CR =Size of original image(bits) Size of compressed image (bits)(1.2) Bits per pixel per band bpppb =Size of compressed image (bits) Lines ×Columns ×Bands (1.3) 1.2.4.1 Lossless versus lossy compression As it was mentioned, compression can be lossless or lossy. Lossless compression allows the source data to be reconstructed perfectly, therefore it is acceptable to be used to compress any kind of scientific data without sacrificing data quality. On the contrary, lossy compression removes some information in order to achieve higher compression ratios. The removed information cannot be recovered when the data are decompressed. The compression ratios for lossless techniques are typically limited to values around 2 or 3, except for data with low information density such as a black sky. Another limitation of lossless compression is the fact that the compression ratios cannot be predicted in advance, what makes it difficult to estimate the amount of data that will be sent, i.e. the necessary downlink capability. Furthermore, it is necessary to introduce error containment Chapter 1. Introduction 15 strategies, since when the data are corrupted, it can no longer be perfectly reconstructed and errors can propagate. Lossy compression allows to find a trade-off between source fidelity and compression ratio and achieves significantly higher compression ratios, although it shares some of the limitations of the lossless techniques, like the possibility of error propagation. There is no upper limit for the ratio that can be achieved with lossy compression techniques, which is established depending on the downlink constraint and the amount of original data measured. Furthermore, progressive compression provides a bridge between lossless and lossy compression methods. It partitions the data into ordered hierarchical segments. Each compressed segment, when combined with the previous ones, allows for the reconstruction of successively higher fidelity versions of the data. The initial version of the reconstruction is very lossy, while the final reconstruction can in principle be lossless or nearly lossless. 1.3 Motivation of research While the resolution of the remote sensors, and consequently the data rates, continue to increase, the available downlink bandwidth is comparatively stable, as has been observed by NASA [27] and European Space Agency (ESA) [28]. The solution offered is to apply data compression, hence payload data processors on-board satellites have to be able to accomplish this task. In particular, Earth Observation missions have the highest performance needs for data processing, data reduction and compression, and future missions and applications will require more powerful on-board processing platforms. 22 Chapter 1. Introduction The validity of the proposed solutions is demonstrated by showing the correctness of the algorithm’s execution on the different platforms, comparing the results with those obtained with golden reference software implementations. Further goals are achieved part of the progress of the work plan. These are detailed next. Different digital hardware design methodologies are utilized and compared, including high-level synthesis solutions. Accelerating algorithms by means of parallelization implies solving bottlenecks and data dependencies which can reduce the parallelization capabilities. The common difficulties which appear when trying to accelerate hyperspectral image compression algorithms are identified. Solutions are given in order to solve data dependencies and allow for more parallelization of the algorithms. When possible, the parallelization strategies are given in a general way, so that they can be applied for other algorithms of the same kind. All these contributions are expected to be useful to reduce the cost and improve the performance of future satellite missions in which hyperspectral on-board data compression will play a critical role. 1.5 Organization of this document The present document is structured in six chapters, including this introductory one, which is dedicated to present the main motivations and goals of this Thesis work. The rest of the chapters are briefly described next. Chapter 1. Introduction 23 Chapter 2: On-board hyperspectral image compression algorithms and hardware implementations This chapter presents a review of the state-of-the-art in the field of algorithms for hyperspectral data compression on-board satellites and their hardware implementations. The main requirements of the algorithms and the on-board hardware are described and the different solutions proposed in the literature are analysed, with the objective of contextualizing this research work and setting quantitative goals. Chapter 3: Implementation of a lossy compression algorithm for hyperspectral images on a GPU This chapter describes the GPU implementation of a lossy compression algorithm for hyperspectral images, showing the acceleration that can be potentially obtained. Chapter 4: Implementation of a lossy compression algorithm for hyperspectral images on an FPGA This chapter shows the implementation on a FPGA of a lossy compression algorithm for hyperspectral images. The resulting experimental results are useful to evaluate how well the algorithm performs on an FPGA and furthermore make it possible to perform a comparison with the GPU implementation of the same algorithm developed in Chapter 3. Chapter 5: Implementation of the CCSDS standard for lossless hyperspectral image compression on a space-qualified FPGA An implementation of the CCSDS standard for lossless hyperspectral image compression is performed on a space-qualified FPGA, which is currently being used in space-applications. Several architectural options are explored for the implementation of the algorithm, selecting the one that provides the best trade-off between complexity, power consumption and flexibility of the 24 Chapter 1. Introduction implementation. Experimental results are provided and a comparison with the other implementations studied in this Thesis as well as implementations of the state-of-the-art is given. Chapter 6: Conclusions Finally, the collection of contributions of this Thesis are summarized and further research works are proposed. Chapter 2 On-board hyperspectral image compression algorithms and hardware implementations This Chapter summarizes the main characteristics of the algorithms for satellite data compression, focusing on those specifically designed for on-board hyperspectral image compression. An overview of the particular requirements of on-board compression is provided, as well as a detailed description of the most relevant algorithms of the state-of-the-art, based on their theoretical basis. Furthermore, the physical implementations and performance figures of several hyperspectral compression algorithms on different hardware technologies, including space-qualified devices, are presented. 25 26 Chapter 2. On-board compression algorithms and hardware 2.1 Outline An efficient compression of hyperspectral images on-board satellites is mandatory in order to save bandwidth and storage space. Hence, it has become a very popular research topic for academia and the space industry. Reducing the data volume in a harsh environment like space, where the computational power is limited, is also a challenge which has been faced with a twofold approach: to propose new algorithms, specifically designed to take advantage of the nature of hyperspectral images; and to present the technologies and strategies to execute the compression in the hardware available on a satellite, minimizing the complexity, and consequently the resource usage and power consumption. In order to make it easier for the space industry to implement on-board compression on their satellites, an effort has been made to develop standard algorithms. In this sense, the Consultative Committee for Space Data Systems (CCSDS), a consortium of the major space agencies in the world, has issued three recommended standards for space data compression: a universal lossless compression solution [24], a lossless to lossy two-dimensional (2D) image compressor [34] and a lossless compression algorithm for multispectral and hyperspectral images [25]. Both, lossless and lossy techniques for hyperspectral image compression can be found in the literature. As any other state-of-the-art compression algorithm, they utilize the redundancies in the image samples to reduce the data volume. Two different approaches are feasible in the case of hyperspectral image compression, either 2D coding, which only takes advantage of the spatial redundancies among neighbouring pixels, or three-dimensional (3D) coding, which also exploits the existing redundancies between bands. Among the 2D approaches it is possible to find algorithms like LOCO-I [35] or 2DCALIC [36]. However, since most remote-sensing images have a large number of spectral components (hundreds of bands in the case of hyperspectral images), taking into consideration the third spectral dimension has been proven Chapter 2. On-board compression algorithms and hardware 27 to increase the compression performance of the algorithms, achieving higher compression ratios. Examples of these are the 3D extensions of the aforementioned algorithms, LCL-3D [37] and 3D-CALIC [38]; and the algorithms presented in Section 2.2. Compression algorithms are inherently computationally demanding, and those designed for hyperspectral images are not an exception. When they are to be executed on a satellite, the specific requirements of on-board processing have to be considered. The available processing power is limited, and most usual data compression technological solutions used on ground cannot be applied to space data systems. Therefore, together with the algorithms, different hardware implementations have been proposed for the compression to be executed on-board a satellite. The most usual approach is to implement a demonstrator of the algorithms in software, which is executed on a general-purpose single-threaded central processing unit (CPU). However, the processors available for on-board usage are not powerful enough to accomplish the compression in an efficient way. Other solutions, more suited for on-board compression, include highly-customized ASICs and implementations on reconfigurable FPGAs. Lately, GPUs have been shown as a promising alternative, although they cannot be used for on-board compression nowadays. A discussion about the current state-of-the-art hardware implementations for hyperspectral image compression is presented in Section 2.3. 2.2 Algorithms for on-board hyperspectral image compression Typically, a hyperspectral compression algorithm consists of a spatial and/or spectral decorrelator, a quantization stage and an entropy coder, which explores the probability of the symbols to assign short codewords to the most 28 Chapter 2. On-board compression algorithms and hardware probable symbols and long codewords to the less probable ones. The decorrelator can be transform-based or prediction-based. In the former approach, a transform like the Discrete Wavelet Transform (DWT), Kahrunen-Lo`eve Transform (KLT), or Principal Component Analysis (PCA) is utilized to decorrelate the data, whereas in the latter, the samples are predicted from neighbouring (in the spectral or spatial direction) samples, and the prediction errors are encoded. Lossless algorithms have been traditionally preferred, to preserve all the information present in the hyperspectral cube for scientific purposes. However, the performance in terms of compression ratio of the lossless techniques is limited, usually showing a data volume reduction from half up to one third the size of the original image, at most. The necessity for obtaining higher compression ratios will become more critical in the near future, as the datarate of the next-generation sensors is expected to increase. Therefore, lossy techniques have been also proposed. Whenever lossy techniques are employed, it is necessary to evaluate the impact of the losses in the reconstructed data. Usually the relationship between the peak signal-to-noise ratio (PSNR) and the compression ratio, commonly named rate-distortion (RD) ratio, is used as a metric, together with the maximum absolute error (MAE) and the mean square-root error (MSRE). The aforementioned metrics are defined in the following equations, where sz,y,x represents a sample located in coordinates (x, y) and band zand ˆsz,y,x is the corresponding reconstructed sample. Ny,Nx and Nz are the number of lines, columns and bands respectively. Compression ratio in bits per pixel per band bpppb =Size of compressed image (bits) Nz ×Ny ×Nz (2.1) Chapter 2. On-board compression algorithms and hardware 29 Maximum absolute error (MAE) MAE = max (|sz,y,x −ˆsz,y,x|) (2.2) Mean-squared error (MSE) MSE = P|sz,y,x −ˆsz,y,x| Nz ×Ny ×Nx (2.3) Peak signal to noise ratio (PSNR) PSNR = 10 log(215 −1)2 MSE (2.4) Nevertheless, it has been observed that a high PSNR does not necessarily yield higher quality in the reconstructed hyperspectral images when they are used in specific applications [39]. Hence, other application-oriented assessments of the impact of the losses have been reported in the literature, where the experiments aim at demonstrating how useful the reconstructed hyperspectral images are at the post-processing stage for particular purposes, e.g. classification, endmember extraction or anomaly detection [40–43]. It is observed that the compression techniques might introduce artifacts which have a little impact in the PSNR but can significantly bias the analysis results of the decompressed image. Lossy compression can also produce a low-pass filtering of the image, which might ease the extraction of the endmembers from the reconstructed data. For instance, in [20] the H.264/AVC video coding standard is utilized to compress hyperspectral images, taking advantage of the fact that both, video and hyperspectral images, constitute 3D data. Despite the high compression ratios achieved, it was demonstrated that a very accurate endmember extraction from the decompressed data was still possible. 30 Chapter 2. On-board compression algorithms and hardware 2.2.1 Requirements and limitations of an on-board hyperspectral image compression algorithm The hyperspectral image compression algorithms need to meet several requirements to be amenable to operate on the hardware available on a satellite, specifically: The complexity of the compression algorithm has to be low. Regardless of the technology where the algorithm is finally implemented, the computational power of the space-qualified hardware used nowadays on satellites is much lower than that of any personal computer or workstation used on ground. Although the radiation-hardened space-qualified components have the same functionality of an equivalent standard processor, they are designed to be insensitive to ionisation and hence are more expensive to design and manufacture. As a result, the available devices do not have the state-of-the-art computational capabilities. Moreover, it is advantageous that the algorithm can be parallelized in order to speed up the compression process for high data-rate sensors; and must use the available resources effectively, possibly not needing an external memory. It is desirable that the algorithm is resilient to errors. Errors can take place during on-board compression because of the effects of radiation, e.g. bit flips in the on-board memory, and while transmitting the data to the ground station. A corrupted packet will prevent the decoding of other packets that depend on it, causing significant error propagation. Traditional compression algorithms cannot recover from a single bit error, causing a wrong decoding of the remainder of the compressed file after the error. Error-resilience aims at limiting error propagation at the cost of losing compression performance, by using error-resilient Chapter 2. On-board compression algorithms and hardware 31 entropy codes or by partitioning the data in units that are coded independently, in such a way that an error in one unit will not prevent from decoding other units. The compression algorithm has to be able to handle raw data. Nevertheless, most of the times the algorithms are only tested on calibrated data. The significance of these results is bounded by the fact that raw data generated on-board are known to have quite different characteristics than calibrated data. For instance, an algorithm for hyperspectral image compression based on look-up tables known as LUT [44] showed remarkable high compression ratios when applied to old calibrated AVIRIS data, because it exploited artificial regularities introduced in the conversion of raw data values to radiance units. These methods did not work that well on raw or newer calibrated data. 2.2.2 Transform-based compression algorithms for hyperspectral images A transform-based compression method applied to 2D data consists of transforming the spatial information to another domain, in such a way that the data are decorrelated. The most popular transforms used in compression are the DWT and the Discrete Cosine Transform (DCT). The transform is followed by the quantization and encoding of the resulting coefficients. Wavelet-based compression techniques are particularly interesting and have shown excellent rate-distortion performance for traditional 2D imagery. Although transform-based methods are mostly popular for lossy compression, reversible transforms allow lossless compression likewise. Examples of popular transform-based compression methods are the JPEG2000 standard [45], and set partitioning methods, such as set partitioning in hierarchical trees (SPIHT) and its 2D and 3D variations (SPIHT-2D, SPIHT-3D, SPECK). 38 Chapter 2. On-board compression algorithms and hardware 2.2.4 Recent research on hyperspectral image compression algorithms Recent work has applied ideas from distributed source coding to construct extremely simple and error-resilient algorithms [67]. Distributed source coding techniques consider a situation in which two or more statistically dependent information sources must be encoded by separate encoders which do not share any information. The theory proves that, under certain conditions, separate coding is optimal, provided that the sources are decoded jointly. When applied to hyperspectral images, the previous band is used for the prediction of the current band. The first band is transmitted uncompressed, while for all others, the prediction parameters are not sent to the decoder. Instead, the decoder reconstructs the pixels by guessing them, and computing a cyclic redundancy check (CRC). Once the CRC matches the one included in the compressed file, the process terminates. This approach provides furthermore error resilience, since an error in the transmitted compressed data does not necessarily yield an erroneous reconstruction. The proposed algorithm shows a competitive complexity when compared with the state-of-the-art, adding the advantage of error resilience features. However, an appropriate trade-off is yet to be found between robustness, complexity and compression performance. The algorithms based on distributed source coding can achieve a compression performance higher than the state-of-the-art 2D prediction algorithms and slightly lower than other 3D prediction-based algorithms, which is the price to be paid for error resilience. Other recent studies are focused on compressed sensing techniques [68], which suggest that a signal, supposed to be sparse, can be perfectly reconstructed from a limited, i.e. fewer than Shannon, number of incoherent measurements. These techniques could indeed simplify the process of hyperspectral image acquisition [69], providing a reduced number of measurements directly produced by the sensor, saving an important amount of resources. Chapter 2. On-board compression algorithms and hardware 39 In fact, preliminary results have demonstrated that the amount of measurements needed to represent a hyperspectral image can be reduced by a factor of up to 10 [70]. Nevertheless, the design of a sensor able to produce these measurements is difficult, and a lot of technological developments are still needed in order to leverage the full potential of this know-how for hyperspectral imaging [71]. 2.2.5 CCSDS Standard algorithms for satellite data compression The importance of an efficient data compression in space missions is further evidenced by the fact that the CCSDS has issued several standards which facilitate for the different the space agencies and industries to exploit the benefits of compression, by making high quality documentation available and helping to establish a broad user community. The first released compression standard for space applications is known as CCSDS 121 [24] and is a universal lossless data compressor consisting of a preprocessor and an entropy coder based on Rice coding [72]. The objective of the preprocessor is to change the statistics of the data by applying a reversible function, hence reducing the entropy. The recommendation does not strictly specify the preprocessing stage, which can be determined by the final user according the specific characteristics of the target data. The subsequent stage consists of a mapper followed by an entropy coder which operates on blocks of Jsamples. It incorporates multiple coding options, based on Golomb power-of-two codes, which are applied concurrently to a J-samples block. Furthermore, it includes a zero-block and a no compression option, as well as a low entropy option known as second extension. The algorithm option that yields the shortest encoded length is selected for transmission. 40 Chapter 2. On-board compression algorithms and hardware The CCSDS 121 features very low complexity, however its performance decreases significantly with the presence of outliers or when the data do not follow any well-defined statistics. With the motivation of overcoming this difficulties, some alternatives to the CCSDS 121 can be found in the literature. In [73] a Fully Adaptive Prediction Error Coder (FAPEC) is presented, together with its software and hardware implementations. FAPEC shows increased compression ratios when compared with CCSDS 121, and a complexity that is amenable for an on-board implementation. More details about the hardware implementation of FAPEC can be found in Section 2.3. The Image Data Compression recommendation [34], CCSDS 122, describes a compression technique which can be used to produce both lossy and lossless compression of 2D satellite images. It consists of a DWT module that performs decorrelation and a bit plane encoder (BPE). The DWT module employs a three-level 2D-DWT, by repeatedly applying a one-dimensional DWT. It is possible to choose between a float DWT or an integer approximation to this transform. The output coefficients are converted to integer values before applying the BPE, which represents each value with a binary word consisting of a single sign bit along with several magnitude bits depending on the bit width of the input image data. The CCSDS 122 is similar to JPEG2000, however it has a a reduced performance which allows for lower complexity and hence low-power hardware implementations. Finally, the CCSDS 123 recommendation [25] defines a payload lossless data compressor that can be applied to multispectral and hyperspectral imagers and sounders. The compressor consists of a predictor and an entropy coder and is based on the FL algorithm [55]. The predictor uses an adaptive linear prediction method to predict the value of each image sample based on the values of nearby samples in a small three dimensional neighbourhood. The residual of the prediction is mapped to an unsigned integer value and encoded with an entropy coder whose parameters are adaptively adjusted to adapt to Chapter 2. On-board compression algorithms and hardware 41 changes in the statistics of the mapped prediction residuals. The standard offers the alternative of using the block-based entropy coder defined in the lossless data compression standard, CCSDS 121. Experimental results in terms of compression ratio for a real hyperspectral and multispectral image corpus show that the CCSDS 123 standard is competitive with other state-ofthe-art algorithms, providing the best trade-off between coding performance and computational complexity [74]. Figure 2.2: CCSDS 123 Recommendation for lossless multiand hyperspectral image compression. 2.3 Physical implementations for on-board compression of hyperspectral images Several algorithms with different complexity and performance features for on-board compression have been proposed in the literature. Nevertheless, it 42 Chapter 2. On-board compression algorithms and hardware is also necessary to provide physical implementations which serve to demonstrate that the algorithms are suited for on-board compression, and that their performance will be maintained when implemented in the on-board hardware. Obtaining a physical implementation is a difficult and time-consuming task, therefore not all the studies about hyperspectral image compression algorithms include them. Many times only the theoretical basis is explained, and the experimental results are obtained with high-level software tools, such as MATLAB. In general, the different implementations found in the literature can be classified based on their technology as follows: Software implementations: - Implementations on general-purpose CPUs. Usually a high-level programming language is used to produce a software which can be executed on any general-purpose CPU. These kinds of implementations are really flexible and can be obtained in a relatively short time at a low cost. However, the low throughput presented by the on-board processors when performing hyperspectral image compression makes them inadequate to perform hyperspectral image compression on-board a satellite, as it will be further explained in Section 3.3. - Implementations on DSPs. These devices are specifically designed to perform signal processing in an efficient way and are therefore an interesting option for hyperspectral image compression, due to the amount of mathematical operations demanded by the algorithms, mainly in those based on transforms. - Implementations on GPUs. These devices have recently become popular for general-purpose computing, employing massively parallel processing to achieve high throughput. Nevertheless, they are Chapter 2. On-board compression algorithms and hardware 43 not space-qualified and their usage on-board satellites is subtle to future technical developments which would reduce their power consumption and make them insensitive to solar radiation. Hardware implementations: - Implementations on ASICs. These highly-customized implementations achieve high throughputs and low power at the cost of relatively high design times, more expensive manufacturing and lack of flexibility. - Implementations on FPGAs. Represent a trade-off between customization and cost, making it possible to obtain high throughputs and low power consumption. The remaining of this Chapter gives more details about the aforementioned technologies and the most significant hardware implementations of several of the algorithms described in Section 2.2. It has to be noted that not all the presented implementations correspond to 3D hyperspectral compression algorithms. FPGA implementations of universal satellite data compression algorithms or 2D image compressors are also presented, since they are considered relevant for the state-of-the-art and they can anyhow be employed for hyperspectral image compression with a reduced compression efficiency, as it was shown in Section 2.2. 2.3.1 On-board hardware technology requirements In general, any hardware implementation, regardless of the technology, is desired to be small in area, have low power consumption and achieve a high throughput. In this sense, the implementation of algorithms for on-board hyperspectral image compression are no exception. However, the specific 44 Chapter 2. On-board compression algorithms and hardware characteristics of the space environment makes it mandatory for the hardware operating on a satellite to meet additional requirements, mainly tolerance to solar radiation and low power consumption, among others. All this difficulties the on-board device fabrication, increases its cost and ultimately reduces its performance. The most relevant requirements of the on-board hardware are summarized next. Tolerance to high energy radiation Hardware operating in space has to be integrated by components and manufactured with materials which can tolerate high energy radiation. The biggest threat to their operation are high energy particles like galactic cosmic rays, solar winds, solar events and radiation belts. High energy particles can cause upsets in complementary metal-oxide semiconductor (CMOS) and field effect transistor (FET) technologies, which are used to manufacture most of nowadays microelectronic devices. A single event upset (SEU) is a change of state caused by ions in a microelectronic device, producing an error in the device output or even a permanent damage which can destroy the device. The relevant effects of the space radiation environment on microelectronics can be divided into two categories: total ionizing dose (TID) effects and single event effects (SEE). The SEEs can be further subdivided into effects that lead to permanent damage (latch-up) and recoverable effects. Radiation effects that accumulate over time are referred to as total ionizing dose (TID). In normal operating conditions, a voltage applied on te gate of a MOSFET transistor creates an electric field, which reaches into the semiconductor below the gate oxide. This electric field causes the formation of a conducting channel between source and drain. When high energy particles impact on the gate oxide of a transistor, they cause ionisation. Electrons get swept out, leaving behind immobile holes, as it is shown in Figure 2.3. These Chapter 2. On-board compression algorithms and hardware 45 positively charged holes in the gate oxide decrease the threshold voltage required for the creation of a conductive channel. If the threshold becomes too low the device is in a permanent “on” state. Space-qualified microelectronic materials are tested for a specific total ionising dose, to ensure they operate as expected in space. Channel On with VG<VT Gate Oxide Drain Gate Source P-Type Silicon VG=0 n+ n+ ++++++ Conductive Inversion Channel (VG>VT) Gate Oxide Drain Gate Source P-Type Silicon +VG n+ n+ (a) (b) Figure 2.3: Radiation effects on a MOSFET transistor.(a) Normal operation. (b) Post irradiation. A latch-up is the occurrence of a path of low resistivity between the voltage supply and ground connection, which can destroy the device. On the other hand, protons can origin nuclear reaction with silicon atoms, what produces short range ionisation and can upset a memory cell. When a heavy ion impacts on material, it loses energy which causes ionisation of the atoms in the proximity of the impact trace. When ionisation occurs in a pn-junction of a microelectronic device the charges created are separated due to the pn-junction’s electric field. This process creates a charge in the electronic device, what can cause an erroneous transient or interact with the charge of a memory cell. If the charge created exceeds a critical value, the state of the memory cell can be upset. When microelectronic devices for space are manufactured with a technology susceptible to SEUs, mitigation techniques have to be employed in order to 46 Chapter 2. On-board compression algorithms and hardware reduce the system error rate. The typical approach is to apply triple modular redundancy (TMR), which uses three replicas of the same circuit and applies a majority voting strategy to select the correct output among the replicas. The direct consequence of utilizing TMR is an increase in the design area. Some commercial off-the-shelf (COTS) solutions have been tried in the past to implement on-board compression algorithms. In particular, in [75] an implementation study of JPEG2000 standard is performed with radiation hardened components, exhibiting disappointing results because of the complexity of the algorithm for its implementation in such a hardware. Subsequently, in the same study, a commercial hardware platform implementing JPEG2000 is used to perform a performance analysis and a study of tolerance to the spatial environment. The circuit revealed a high sensitivity to radiation, demonstrated by the fact that not a single image was compressed successfully under heavy ions beams. Although the COTS solution is really efficient to implement JPEG2000, it is really inappropriate for on-board usage. Examples like these illustrate the need for hardware specifically designed to operate on-board a satellite. Power efficiency Besides the device being radiation-tolerant, it has to be considered that the available power on-board a satellite is limited and therefore the constraints in terms of power consumption for on-board technologies are much more restrictive than for commercial applications. High throughput In the specific case of hyperspectral image compression, taking into account the amount of data to deal with, a high throughput is desired, especially when real-time compression is wanted. Although compression is expected to become a necessity in the future space missions, there is still a lack of viable on-board platforms to perform significant image processing and compression. Chapter 2. On-board compression algorithms and hardware 47 2.3.2 Software implementations A hyperspectral image compression algorithm can be described by a set of instructions suitable to be executed in a computer’s processor. Among the existing software implementations, it is possible to find implementations on general-purpose processors, digital signal processors (DSPs), and graphics processing units (GPUs). The main contributions that can be found in the literature are summarized in the following sections. Moreover the advantages and disadvantages of each approach are explained, from the point of view of the on-board hyperspectral image compression needs. 2.3.2.1 Implementations on general-purpose CPUs Most of the algorithms proposed for hyperspectral image compression have been implemented in software for general-purpose processors. The programming languages utilized vary, being C/C++ and Java the most popular. These software implementations can be compiled and executed on any processor, and most of them are open source. Their purpose is usually to demonstrate how efficient the algorithms are in terms of compression ratio. They are also sometimes utilized to address the complexity of the algorithms, by comparing for instance the execution times of two different algorithms. Although it can be in general inferred that the faster the execution of the algorithm, the lower its complexity, the execution times are also dependent on how the software was programmed -e.g. to make the source code understandable or to maximize the throughputor the compiler options. Hence, in order to address the complexity of an algorithm other metrics besides the software execution times shall be used, for instance the number of operations, the bit width of the different variables, the precision of the operations or if they are integer or floating-point operations. 54 Chapter 2. On-board compression algorithms and hardware determined. Afterwards, the behaviour of the circuit is described. This description can be performed at different levels of abstraction. The most commonly utilized is register-transfer level (RTL). At RTL, a synchronous design is modelled in terms of logical operations and the data path between hardware registers. Hardware description languages like Verilog or VHDL are utilized in order to create a high-level representation of a circuit, from which a lower level representation and actual wiring (physical design) can be derived, performing the logic synthesis, placement and routing. Finally, the design is prepared for manufacturing depending on the final implementation technology. The RTL design flow is rather long, and it involves careful planning and a designer experienced with hardware description languages. Verification is necessary in order to match the behaviour of the design with the specification at various stages of the design flow. The process can be shortened with the help of high-level synthesis tools, which are capable of porting high-level source codes written in programming languages like C/C++ or Matlab to RTL. Examples of these are CatapultC, Impulse C or C-to-Silicon [83]. These tools are useful for rapid prototyping, and present several advantages, e.g. the amount of code to be written by the designers is highly reduced and the tools open up opportunities for extensive design space exploration. However, if a highly optimized implementation is desired, it is mandatory to write the RTL description from scratch. 2.3.3.2 Implementations on ASICs ASICs are tailored to be optimum for a particular application and therefore can achieve high throughputs and low power when specifically designed for on-board hyperspectral image compression. Most of the 2D compression algorithms which are currently executed on satellites are implemented on ASICs [84]. Chapter 2. On-board compression algorithms and hardware 55 An example of an ASIC implementation for space is the CWICOM chip, see Figure 2.6, developed by Astrium in the frame of an ESA contract. CWICOM is a high speed 2D image compression ASIC which implements the CCSDS 122 standard for image data compression [34]. The CCSDS 122, as it was explained in Section 2.2.5, consists of a DWT followed by a BPE and can perform lossless as well as lossy compression. Figure 2.6: CWICOM compression ASIC. CWICOM features high data rate of 60 Msamples/sec, what results in 960 Mbps if the sample values are represented with 16 bits. The power consumption is <100 mW/Msample/sec. This means that for the maximum throughput of 60 Msamples/sec the consumption would be approximately 3 W. It supports a TID of 100 KRad and is tolerant to SEU thanks to its internal EDAC. One of the main challenges present during the development of a hardware implementation of the CCSDS 122 algorithm is the high amount of memory which is needed to store the DWT coefficients during the processing. An external memory can address this problem, but would decrease the performance. Hence, the CWICOM includes a high amount internal memory cells, and a very efficient internal embedded memory organization, making it possible to compress the images without the need of an external memory. While ASICs are power and area efficient, they also present some weaknesses. In many occasions, space applications require more flexibility and 56 Chapter 2. On-board compression algorithms and hardware scalability for post-launch modifications and repair. ASICs, once manufactured, cannot be configured to efficiently match subsequent mission needs and requirements. 2.3.3.3 Implementations on FPGAs Alternatively to ASICs, FPGAs offer solutions at a lower cost and with increased flexibility. They consist of arrays of blocks of generic logic cells which can be interconnected in a general way. Both the logic blocks and the interconnection structure are programmable by means of switches which can be set as open or short circuit. FPGAs are manufactured by a number of companies like Microsemi, Altera or Xilinx, among others. The basic scheme and elements of an FPGA are shown in Figure 2.7. I/O Pad Routing Line Logic Resource RAM RAM RAM RAM Combinational Sequential Figure 2.7: Basic elements of an FPGA. FPGAs present several advantages for on board hyperspectral image compression. They are able to apply parallel processing to increase throughput and provide some flexibility at the same time. Moreover, FPGAs offer the Chapter 2. On-board compression algorithms and hardware 57 possibility of adapting the designs to successive upgrades of compression algorithms and electronic components over a long term. Moreover, FPGAs provide scalability and data integrity features. The switches which allow the FPGA programming can be constructed in several ways, including: pass-transistors constrolled by RAM, anti-fuses, EPROM transistors or EEPROM transistors (flash-based). Both static random access memory (SRAM)-based and anti-fuse FPGA are used nowadays in satellites. Traditionally, anti-fuse technology has been preferred by the space industry, due to its increased robustness against radiation. An anti-fuse normally resides in a high-impedance stage, but can be “fused” into a low-impedance state when programmed by a high voltage. An example of radiation-tolerant anti-fuse FPGAs is the RTAX family from Microsemi. However, anti-fuse FPGAs present an important disadvantage: they are one-time programmable. On the other hand, in SRAM-based FPGAs the programmable connections are made using pass-transistors, transmission gates or multiplexers that are all controlled by SRAM cells. The major advantage of this technology is that it provides an FPGA that can be reprogrammed many times and very quickly, and can be produced using standard CMOS technology. However, one of their disadvantages it that they have to be programmed again every time the system is powered up, due to the volatile nature of SRAM. SRAMbased FPGAs are becoming increasingly important for new space missions. The lifetimes of the satellites expand far beyond 10 years, which is much longer than the validity of telecom standards. Hence, reprogrammability becomes an important requirement. 58 Chapter 2. On-board compression algorithms and hardware SRAM-based FPGA are susceptible to SEUs, which are induced by highenergy particles in the harsh environment of space. This problem is addressed through the use of radiation-tolerant and radiation-hardened technologies as well as SEU mitigation techniques. Unlike ASICs, where only memory elements have to be protected against SEU, SRAM-based FPGAs must implement full triplication of all the design elements due to the high sensitivity of its configuration memory to radiation effects. A single error affecting a TMR circuit is masked and tolerable, however, when multiple independent SEUs hit a TMR circuit, it might occur that the majority voter in the TMR scheme votes for the wrong answer. To mitigate his effect, the configuration bit-stream has to be written periodically back to the FPGA, whith a scrub cycle that must be selected according to the frequency of SEUs occurrence. During this reconfiguration time, which can be of the order of milliseconds, the device is offline, what can be intolerant for systems with hard real-time constraints. The FPGA manufacturer Xilinx has developed several radiation-hardened models as part of the Virtex IV and Virtex 5 families. The radiationhardened Xilinx FPGA XQR40662XL on-board FedSat, is the first demonstration of hardware reconfiguration in space [85]. Finally, flash-based FPGAs are similar to SRAM-based FPGAs, but they present the advantage that the configuration memory is non-volatile and therefore they can be programmed off-line. The downside is that flashbased FPGAs tend to have higher static power consumption and longer programming times. An example of radiation-tolerant flash-based FPGAs are the ProASIC3 family, manufactured by Microsemi. This model, however, is not used in any current space mission. The most relevant contributions to the field of FPGA implementations of satellite data compression algorithms are summarized next. Chapter 2. On-board compression algorithms and hardware 59 The FAPEC universal lossless data compressor was implemented on a radiationhardened RTAX anti-fuse FPGA, from a RTL description in VHDL language [73]. Several considerations are taken into account in order to lower the complexity of FAPEC in such a way that its FPGA implementation is more efficient. For instance, floating-point operations are avoided and the data block size is reduced. These changes are reported to have little or no effect in the performance of the algorithm. The designed compressor is finally implemented in a ProASIC3L developement board which includes a M1A3P1000L FPGA, whose number of equivalent gates is equal to that of the space-qualified anti-fuse RTAX1000S. The implementation shows a critical path of 18.32 ns (maximum frequency 55 MHz), with a throughput of 2 Msamples/sec (32 Mbps if the input samples are 16 bits wide). The hardware occupancy is 12% of combinational logic and 15% of sequential logic and the power consumption is estimated at 35 mW. Another universal compression algorithm and the corresponding hardware architecture are presented in [86]. The algorithm performs context-based statistical lossless compression of multiple types of data. It takes advantage of FPGAs that support partial and dynamic reconfiguration, which consists of changing a portion of the reconfigurable hardware while the rest is still operating. Their method is based on a dynamically reconfigurable modelling stage followed by statistically configured probability estimation and an arithmetic coding. The dynamic modelling is specialized to each data type and uses a combination of context modelling, predictive coding and motion estimation depending on the data type being processed. The throughput performance of the proposed system is 100 Mbps on a Xilinx Virtex4 SX35 FPGA. Partial dynamic reconfiguration has become a popular research topic in the space industry, since it can enhance space applications with run-time adaptive functionality, enabling mission specific adaptability. However, technical developments are still needed for this technology to become reliable and fault tolerant for its usage on-board satellites [87]. 60 Chapter 2. On-board compression algorithms and hardware Table 2.1: Results of implementing CCSDS 122 on an RTAX2000S FPGA Resource Proba-V EnMap Combinational C-Cells 8455 (39%) 9772 (45%) Sequential R-Cells 7125 (66%) 7620 (71%) Total Cells 15580 (48%) 17392 (54%) Block RAMs 54 of 64 58 of 64 An FPGA implementation of the CCSDS 122 standard for 2D satellite image compression can be found in [88]. The study presents the designed hardware architecture and its implementation on a space-qualified FPGAs, specifically on an anti-fuse RTAX2000S. However, the ProASIC3E flashbased FPGA is utilized for prototyping. As it was explained in Section 2.2.5, the CCSDS 122 standard algorithm consists basically of a Discrete Wavelet Transform (DWT) plus a bit plane encoder. The intermediate coefficients of the DWT have to be stored and rearranged during the compression process, what requires more memory than available in the RTAX2000S FPGA. Hence, an external SDRAM memory is used, with an appropriate memory organization to reduce the access overhead. The hardware occupancy and frequency of the implementation depends on the selected configuration parameters of the CCSDS 122 algorithm. The results of two different configurations for two specific space missions, namely ESA Proba-V and the German EnMap are shown in Table 5.23. The maximum frequency is 64 MHz for the Proba-V mission and 50 MHz for the EnMap mission. The throughput is 90 Mbps for Proba-V and 130 Mbps for EnMap. As it was already mentioned, transform-based methods for hyperspectral image compression require a fair amount of processing memory, are computationally intensive, and not amenable to parallelization due to the sequential nature of the transforms. Hence they are not popular candidates for an FPGA implementation. However, efforts have been made to accelerate Chapter 2. On-board compression algorithms and hardware 61 transforms in order to be able to implement them efficiently on FPGAs, as shown in [89], where the KLT is investigated. A comprehensive analysis of the computations needed to calculate the KLT is performed to inspect the feasibility of different acceleration techniques, such as parallelism. The proposed designs are implemented in a FPGA-based System-on-Chip (SoC), which incorporates a flash-based FPGA, a 32-bit ARM Cortex M-3 microcontroller subsystem and an analogue computing engine. The hyperspectral images are divided in clusters of 32 bands to reduce the processing time and memory requirements. The results show an improvement of more than 54 % of the execution time for the proposed architecture and a power consumption of 225 mW for a cluster of 32 bands. These results demonstrate that it is possible to reduce the complexity of the KLT transform to target space applications. The implementation of the ICER-3D algorithm [90] is another interesting case in which a transform-based hyperspectral compression method is implemented on an FPGA. Specifically, it was implemented on a Xilinx Virtex2 Pro (XC2VCP70) with an embedded PowerPC processor and on-chip bus architecture. The implementation features efficient utilization of offchip memory through internal buffering to minimize intensive input/output operations. The results display a maximum frequency of 50 MHz and a throughput of 4.5 Msamples/sec. The power consumption of the prototype is less than 6.5 W and the hardware occupancy is around 60%. Among the FPGA implementations of prediction-based algorithms for hyperspectral image compression, the one of the FL algorithm is particularly remarkable [55]. It has to be noted that the FL algorithm was selected by the CCSDS for standardization and it represents, with little modifications, the base of the current CCSDS 123 standard for lossless multiand hyperspectral image compression. The FL algorithm was implemented on a Xilinx Virtex4 LX160 SRAM-based FPGA [32]. The implementation has a critical path of 62 Chapter 2. On-board compression algorithms and hardware Table 2.2: Virtex IV LX160 device utilization of the FL algorithm Device Virtex IV LX160 Resource Used Slice 67584 (5%) FIFO/RAMB16 9 (3%) DSP48 6 (6%) Table 2.3: Virtex4 and Virtex2 device utilization of a lossless and nearlossless hyperspectral compression algorithm Device XQR4VLX200 XQ2V3000 Used LUT 10306 (5%) 10248 (35%) RAMB16S 21 of 336 (5%) 21 of 96 (22%) Mult18×18 - 9 of 96 (9%) DSP48 9 of 96 (9%) - 29.5 nsec (maximum frequency 33 MHz) and compresses one sample every clock cycle, which results in a throughput of 33 MSamples/sec. It processes the samples in band interleaved by pixel (BIP) order. It is 33 times faster than a software implementation running on a Pentium IV machine. The implementation has a rather low hardware occupancy, as shown in Table 2.2, and its power consumption is estimated at 1.27 W. Another example of a prediction-based algorithm is the lossless to nearlossless hyperspectral compression algorithm presented in [56]. The algorithm implemented on a Xilinx Virtex2 V3000 and on a Xilinx Virtex4 LX200 FPGA, producing the hardware occupancy results presented in Table 2.3, with a maximum frequency of 81 MHz and a throughput of 70 Msamples/sec for the Virtex4 FPGA and 79 MHz and 69 Msamples/sec for the Virtex2. Table 2.4 summarizes the most relevant figures in terms of throughput, frequency and power of some of the FPGA implementations of satellite data compression algorithms presented in this Section. Since the algorithms were Chapter 2. On-board compression algorithms and hardware 63 implemented on different technologies, the hardware occupancy data could mislead the reader to imprecise comparisons and therefore are not included in Table 2.4. Nevertheless, the on-board compression algorithms implemented on FPGAs that can be found in the literature have a complexity which is low enough to allow their implementation on most of the currently available space-qualified FPGAs. 70 Chapter 3. Lossy hyperspectral compression on GPU s0,m,n s0,m,n-1 s0,m-1,n sz,m,n sz-1,m,n band 0 band z band z-1 (a) (b) z m n Figure 3.2: Prediction neighbourhood of the LCE algorithm: a) first band b) all other bands. ˆs0,m,n = (˜s0,m−1,n + ˜s0,m,n−1)1 (3.1) In the previous equation, ˆsdenotes the predictor, ˜sthe decoded value and stands for right shift. For all the other bands, the samples sz,m,n are predicted from the decoded samples ˜sz−1,m,n in the co-located block in the previous band. A least-square estimator (LMS) is computed over the block as α=αN/αD(3.2) where αN=X m,n [(˜sz−1,m,n −µz−1)(sz,m,n −µz)] (3.3) and αD=X m,n [(˜sz−1,m,n −µz−1)(˜sz−1,m,n −µz−1)] (3.4) Chapter 3. Lossy hyperspectral compression on GPU 71 where µzand µz−1are the average values of the co-located decoded blocks in bands zand z−1. Quantized versions of αand µz, denoted α0and µ0 z, are generated using a scalar quantizer. Finally, the predicted values are computed for all samples in a block as: ˆsz,m,n =µ0 z+α0(˜sz−1,m,n −µz−1) (3.5) and the prediction error is calculated as: ez,m,n =sz,m,n −ˆsz,m,n (3.6) 3.2.2 Rate-distortion optimization Before proceeding with the quantization of the prediction error samples, the algorithm checks if the prediction is so close to the actual pixel values that it makes sense to skip the encoding of the prediction error samples and, instead, rise a one-bit-flag indicating that the current block contains all-zero prediction error samples (this is denoted as zero block condition). To make this decision, the energy of the predictor error is computed: D0=1 N×NX m,n e2 z,m,n (3.7) If D0< DT, with DTa chosen threshold, then the zero block condition is triggered and the decoded values are calculated as: ˜sz,m,n = ˆsz,m,n (3.8) 72 Chapter 3. Lossy hyperspectral compression on GPU 3.2.3 Quantization and mapping For the non-zero blocks, the prediction error samples are quantized to integer values e0 z,m,n, and dequantized to reconstructed values e00 z,m,n. The decoded pixel values can be then obtained by adding the dequantized prediction error to the predicted value: ˜sz,m,n = ˆsz,m,n +e00 z,m,n (3.9) For the first band, the quantization process is performed pixel by pixel using a scalar uniform quantizer. For the other bands, it is possible to choose between a scalar uniform quantizer and the uniform-threshold quantizer (UTQ) described in [91]. The UTQ quantizer exhibits superior reconstruction performance, however, this quantizer, as described in the LCE algorithm, cannot be implemented with integer arithmetic only and involves a non-linear formula. This has to be taken into account when implementing the algorithm in hardware. As an alternative, the user can opt for the scalar quantizer, which utilizes only integer arithmetic. Finally, the quantized prediction errors are mapped to non-negative values and entropy encoded, with one exception: in the first band, the first sample of each block is neither mapped nor encoded. It is written in binary format using 16 bits. 3.2.4 Entropy coding The 16 ×16 residuals of a block are encoded in raster order using a Golomb code whose parameter is constrained to a power of two, except for the first sample of each block, which is encoded using an exponential Golomb code of order zero. Chapter 3. Lossy hyperspectral compression on GPU 73 Exponential Golomb (exp-Golomb) codes were first proposed in [57]. The codewords consist of a prefix and a suffix. The prefix part of the exp-Golomb code of parameter kfor a nonnegative integer τconsists of a unary code corresponding to the value u=log2(τ 2k+ 1). The suffix is computed as the binary representation of τ+2k(1 −2u) using k+lsignificant bits. In the particular case of the LCE, the kvalue for the exp-Golomb code is always fixed to 0. The rest of the samples are encoded using a Golomb code. The parameter kjfor the j-th sample of the block is computed from a running count Rjof the sum of the magnitude of the last 32 unmapped residuals of the block; for samples with index less than 32 only the available values are used. Similarly to a Huffman code, Golomb codes are unable to produce codewords shorter than 1 bit. This means that, without the rate-distortion (RD) optimization stage described in Section 3.2.2, it would be impossible to obtain bit-rates lower than 1 bpp. However, with the proposed scheme, the minimum rate for each block is still 1 bpp, except for the blocks that are skipped, what makes it possible for the average rate to go as low as desired. 3.2.5 File format The compressed file is a concatenation of coded blocks which are read spatially in raster order, and each block is coded with all bands. The first band is conformed by the quantized first sample represented with 16 bits followed by the codewords of the remaining mapped prediction residuals of the block. For all other bands, since the predictor is not causal, it is necessary to write the parameters α0and µzfor each block. 10 bits are used to store α0and 16 bits are reserved for µz. These parameters are followed by the zero block decision bit. For the non-zero blocks, the exp-Golomb and Golomb encoded 74 Chapter 3. Lossy hyperspectral compression on GPU codewords are appended. Figure 3.3 illustrates the format of the compressed file. block 0 [ band 0 -> first sample in binary + encoded prediction errors band 1 -> + + zero_block flag + encoded prediction errors ... band Nz] block 1 ... block NB-1 block 0 [ band 0 -> first sample in binary + encoded prediction errors band 1 -> + + zero_block flag + encoded prediction errors ... band Nz] block 1 ... block NB-1 Figure 3.3: File format of the LCE compressed data. 3.2.6 LCE compression efficiency The compression efficiency of the LCE algorithm was compared in terms of RD in [33], showing the relationship between the achieved compression bitrate and the PSNR of the reconstructed images. The results demonstrate that the LCE algorithm shows equal or better rate-distortion (RD) than JPEG2000 with a spectral wavelet transform, but with a significantly lower complexity and memory requirements. The RD curves obtained for an hyperspectral image taken by AVIRIS are depicted in Figure 3.4. The Figure is shown exactly as it appeared in [33]. 3.3 Software implementation of the LCE algorithm The LCE algorithm was originally implemented in software using C programming language, to be executed on a single threaded CPU. It operates independently on every N×Nblock with all its bands. For the specific case of this study, we consider the block size N= 16, which as stated in Chapter 3. Lossy hyperspectral compression on GPU 75 Figure 3.4: Rate-distortion curves for AVIRIS when compressed with LCE and other algorithms of the state-of-the-art [33] [33] typically optimizes the algorithm performance. Each block in the image is identified with two coordinates in the vertical and horizontal spatial directions, namely vb and hb.Bz,vb,hb represents a 16 ×16 pixel block in coordinates (vb, hb) and band z. A single sample in Bz,vb,hb is symbolized as sz,m,n. Additionally, Nvb and Nhb are the total number of 16 ×16 blocks that can be found in the image, in vertical and horizontal spatial dimensions. The LCE compressor consists essentially of a chained loop, which iterates to cover all horizontal and vertical blocks in the image, and all bands in a block. The innermost loop reads a 16 ×16 block Bz,vb,hb and performs the different stages of the LCE algorithm presented in Section 3.2 for every single pixel sz,m,n, namely the prediction, rate-distortion optimization, quantization, entropy coding and bit packing to create the compressed file. The pseudo-code 76 Chapter 3. Lossy hyperspectral compression on GPU 2-D predictor Spectral predictor Rate-dist optimization Quantization and mapping Golomb entropy coding Bit packing Band++ 1st band? yes no no yes Zero block? Figure 3.5: Flowchart of the LCE algorithm showing the chained loop and the different compression stages is shown in Figure 3.6. 3.3.1 Generation of the compressed file The result of the previously described chained loop is a bit stream which represents the compressed image. The codewords generated at the entropy coding stage and the parameters α0and µ0 zare written by the software in a single file, in the same order they are obtained. Since Golomb codes produce Chapter 3. Lossy hyperspectral compression on GPU 77 Additional for loops are needed to cover all the pixels in a block Figure 3.6: Pseudo-code of the main function of the LCE algorithm implementation in C language. codewords of variable length, they are buffered in a bit-by-bit fashion in a 32-bit variable. When the buffer is full, it is written to the compressed file, generating the bit stream. 3.3.2 Configuration parameters The software implementation allows the user to configure the LCE algorithm by selecting several parameters to set the functionality mode (baseline or advanced) and the quality of the resulting compressed image. The baseline algorithm performs compression using a scalar quantizer, fully implemented with integer arithmetic; and an advanced algorithm replaces the scalar quantizer with the uniform-threshold quantizer (UTQ) [91], which uses floating-point arithmetic. The parameter UTQ can be set by the user to select the functionality mode. Setting UTQ = 1 enables the UTQ quantizer and setting UTQ = 0 uses the baseline functionality mode. 78 Chapter 3. Lossy hyperspectral compression on GPU The quality of the resulting compressed image can be likewise set by the user by assigning values to a parameter named delta, making it possible to find a trade-off between quality and bit rate. Delta sets the quantization step size of the quantizer, therefore, increasing delta yields to higher compression ratios, but lower quality of the reconstructed image. Delta has to be an integer greater than 1, with delta = 1 providing the maximum possible quality. Taking into account the specifications of the LCE algorithm, we can conclude that it can benefit from parallel execution schemes, due to the amount of operations that can be performed independently. As it was explained before, the algorithm compresses independent blocks of data, and having a closer look at the stages of the algorithm, we observe that there are operations that can be performed on every sample independently. Given the amount of data that can be processed in parallel, we first explore how much the algorithm can be accelerated when executed in a GPU. In the following, we present a brief description of the GPUs architecture and programming environment. Afterwards, we show the main strategies utilized to optimize the LCE implementation in order to achieve the maximum possible speedup with the GPU. Finally, we show the achieved acceleration with experimental results running on real hyperspectral images. 3.4 GPU architecture and NVidia CUDA GPUs are structured as a set of multiprocessors, each composed of a simple processing element working in single-instruction multiple data mode. They offer the possibility to dramatically increase the computation speed in applications where a huge amount of data can be processed in parallel. At the time this document is written, the biggest GPU manufacturers are Intel, Chapter 3. Lossy hyperspectral compression on GPU 79 Nvidia and AMD. GPUs were initially designed to accelerate graphic computations, but recently they have become popular also for general-purpose computing, i.e. they can be programmed to compute any computable value. However, they are very restrictive in operations and programming, and are only effective for problems which can be solved using stream processing - processors that can operate in parallel by running one kernel on many records in a stream at once; and the hardware can only be used in certain ways. The currently dominant open general-purpose GPU programming language is OpenCL [92], but other proprietary frameworks exist. The GPU utilized for this Thesis work is the Tesla C2075 from Nvidia (see 3.4.3), and consequently, it was decided to program it using Nvidia’s platform, known as CUDA. The Computer Unified Device Architecture (CUDA) [93] [94] [95] is a parallel computer architecture developed by Nvidia, which provides a scalable programming model and a software environment for parallel computing. It is an extension of C language, which offers an Application Programming Interface (API), allowing programmers to use a GPU as a massively multi-threaded general purpose co-processor. The APIs make it possible to manage devices and memories, but hide the real hardware from the developers, making it unnecessary for the programmers to explicitly manage threads. The CUDA abstractions guide the programmer to partition the problem into coarse subproblems that can be solved independently in parallel and then into finer pieces that can be solved cooperatively in parallel. 3.4.1 CUDA abstractions CUDA names the CPU and its memory the host and the GPU and its memory, the device. The code executed in parallel using the GPU is typically called a kernel. A kernel is launched from the host and can have parameters, like any other C language functions (see Figure 3.7). 86 Chapter 3. Lossy hyperspectral compression on GPU or spectral dimension, achieving high performance at the same time. Hence, the number of CUDA threads and blocks called in each kernel is established in such a way that the CUDA application is suitable to compress any hyperspectral image regardless its spatial or spectral size in short compression times. Nevertheless, the kernel dimensionality might be further optimized for a specific image, taking into consideration its size and adjusting the number of threads and blocks accordingly. 3.5.1 Allocation of the image data in the GPU The necessary image data has to be sent and stored in the GPU dedicated memory in order to make it possible to perform operations on them. It is decided to copy the whole hyperspectral cube in the GPU before executing any of the kernels, in order to minimize the number of transactions between the CPU and GPU. Once the necessary data are stored in the GPU, an efficient use of the different memory spaces is necessary to hide latency. When data are copied from the CPU to the GPU, they are initially stored in global memory. However, CUDA local or shared memory spaces are used whenever it is convenient in order to accelerate memory accesses. The hyperspectral data are read in the same format which was used in the C implementation, i.e. band-sequential (BSQ). Once allocated in the GPU, the different 16 ×16 pixel blocks are read and copied to the shared memory to hide the global memory latency, since the original pixels are accessed repeatedly and shall be visible for all threads in a block. 3.5.2 Prediction, quantization and mapping The main difficulty found when parallelizing these stages of the LCE algorithm is that it is impossible to avoid having a sequential loop to cover Chapter 3. Lossy hyperspectral compression on GPU 87 all bands, because the processing in band zcannot be started until the processing of samples in band z−1 is finished. Acceleration by means of parallelization is obtained by taking advantage of the fact that each 16 ×16 image block with all its bands can be processed independently. The parallelization strategy, as shown in Figure 3.11, is based on mapping the problem to the CUDA abstractions in such a way that a CUDA block is launched by the kernel for each spatial 16 ×16 block Bvb,hb with all its bands. Each CUDA block is set up to have 256 (16 ×16) CUDA threads, each responsible for performing operations on a sample sz,m,n in parallel. Therefore, the number of CUDA blocks which are launched can be calculated for a specific hyperspectral cube as: CUDA blocks =Ny ×Nx 16 ×16 (3.10) where Ny is the number of lines and Nx represents the number of columns in a specific band of a hyperspectral image. Nz 16x16 block CUDA block CUDA grid 16x16 block Figure 3.11: CUDA abstractions for the parallel execution of the LCE prediction, quantization and mapping stages. 88 Chapter 3. Lossy hyperspectral compression on GPU The hyperspectral image samples are stored in the CUDA global memory. As the samples are going to be accessed many times by the CUDA threads, it is decided to copy them to the CUDA local memory. Specifically, every iteration of the loop that covers all bands, a spatial 16 ×16 block is copied to the CUDA local memory, indexing it with the CUDA blocks and threads identifiers. This way, threads have fast access to the sample they are entitled to process. Local memory can only be accessed by a single thread, and it is invisible for the other threads in a CUDA block. Nevertheless, not all the operations on a sample sz,m,n can be performed independently, as it is the case of computing parameters αand µ, for which a summation has to be performed. However, these operations can be accelerated by making threads cooperate and share information by means of CUDA shared memory. Specifically, in order to reduce the number of iterations of the loops when computing summations, the so-called parallel reduction algorithm is employed [97]. The strategy presented in [98] is also utilized to reduce the number of iterations when prefix-sums have to be performed. Once the prediction errors are obtained, they are copied again to the CUDA global memory, to make them accessible for the next stages of the LCE. 3.5.3 Entropy coding Unlike the prediction stage previously described, the entropy coding operations can be performed on every 16 ×16 block of a specific band in parallel, without any information from neighbouring bands, what makes it possible to process more data in parallel. A kernel is designed to perform the entropy coding of the mapped prediction residuals, which processes each 16 ×16 block of prediction residuals in parallel, as shown in Figure 3.12. For optimization purposes, it is established Chapter 3. Lossy hyperspectral compression on GPU 89 that each kernel launches the maximum number of possible threads allowed by the Tesla C2075 GPUs, which is 1024. Therefore, the number of CUDA blocks to be launched can be calculated as: CUDA blocks =Nz ×Ny ×Nx 16 ×16 (3.11) 4 spatial blocks of 16x16 samples CUDA block CUDA grid Figure 3.12: CUDA abstractions for the parallel execution of the LCE entropy coding stage. It has to be considered that the maximum number of CUDA blocks is limited to 65535. If the number of CUDA blocks calculated in the previous equation is higher than this maximum, then the kernel has to be called more than once, affecting the performance negatively. Setting the number of threads to the maximum is likewise a way to minimize the number of necessary CUDA blocks, and therefore reduce the impact of having to invoke the kernel repeatedly. Several facts are considered when designing the GPU implementation of the entropy coding stage: 90 Chapter 3. Lossy hyperspectral compression on GPU The codewords shall be computed in parallel for each mapped prediction error sample. It is not possible to directly write on a file from the GPU, consequently, the codewords have to be saved to variables in the CUDA memory spaces. The codewords have variable lengths, which can be greater than 32 bits. The strategy followed to generate the codewords and pack them to a bitstream must be different to the one followed in the CPU implementation, which is based on the sequential ordered generation of the codewords, as it was explained in Section 3.3. For the GPU implementation, the approach is to pre-process the Golomb parameters of all mapped prediction error samples in a 16 ×16 block, and afterwards compute the codewords for every sample in parallel. Computing the codewords A strategy is designed to compute the Golomb parameter of every j-th sample of the block, kj, in parallel. This parameter is computed from a running count Rjof the sum of the magnitude of the last 32 unmapped prediction errors of the block, ej; for samples with index less than 32, only the available values are used. This implies that the Golomb parameter of a specific prediction error in a block depends on the accumulated sum of the previous unmapped prediction errors. The running count is calculated for every sample as: Rj=Rj−1− |ej−33|+|ej−1|(3.12) Chapter 3. Lossy hyperspectral compression on GPU 91 In order to be able to obtain Rjin parallel, the prefix-sum of all the unmapped values is computed as: Ej=     0,if j= 0 254 P j=1 ej,if j > 0(3.13) This is implemented in CUDA following the scheme proposed in [98], as it is explained in the following. We note that this strategy requires the cooperation of threads, therefore the necessary data has to be copied to the CUDA shared memory. The array of unmapped prediction errors in a block of 16×16 is represented as: [e0, e1, ...e255] The prefix-sum gives the result: E=                         E0 E1 E2 ... Ej ... E255                         =                         0 |e0| |e0|+|e1| ... |e0|+|e1|+... +|ej−1| ... |e0|+|e1|+... +|e254|                         (3.14) After the prefix-sum, Rjis calculated by subtracting element Ej−33 to each Ejwith j > 32 in parallel, resulting in the desired Rj: 92 Chapter 3. Lossy hyperspectral compression on GPU Rj=   Ej,if j≤32 Ej−Ej−33,if j > 32 (3.15) Once kjis known for every prediction error sample, the codewords can be created in parallel by the CUDA threads. Each thread computes and saves a codeword in its local memory. The codewords are of variable length and each of them is saved in a 32-bits unsigned integer variable. The Golomb codes can produce codewords of any length, i.e. it might happen that a codeword is longer than 32-bits. For those specific unusual cases, an auxiliary array is created to save the codeword in more than one 32-bits variable. Utilizing this auxiliary array can cause a performance penalty. As part of a future work, this problem can be solved by setting a maximum codeword size so that, in case the codeword is longer than the maximum, the data are not encoded but saved in binary. The size in bits of the codewords is also saved by every thread, in order to be able to create the encoded buffer, as it is explained in the following. Generating a compressed bit stream for each 16 ×16 block The encoded prediction residuals have to be saved in raster order to produce the final encoded 16 ×16 block, which contains the ordered sequence of codewords, without leaving any bits unused between them. The strategy presented in [99] is followed to write every codeword of a 16 ×16 block in a single output buffer, as shown in Figure 3.13. First, the final position of a codeword in the output buffer is calculated. This final position is given by two coordinates: the word position where the codeword is saved, word position; and the bit position in that word where the codeword starts, starting bit of every codeword. These coordinates can be calculated in parallel for every codeword. The number of bits taken by each codeword is known as a result of the entropy coding stage. In the following, Chapter 3. Lossy hyperspectral compression on GPU 93 Data structure of codewords after Golomb Coding codeword 0 Shifted codewords word_position = 0 codeword 1 word_position = 1 codeword 2 word_position = 2 codeword 3 word_position = 3 codeword 255 word_position = 255 word 0 word 1 word pword 2 Data structure of compressed buffer for every 16x16 block OR Figure 3.13: Parallel generation of a compressed 16 ×16 block. we denote qjthe number of bits taken by each codeword j. A prefix-sum of these data yields the bit position of a codeword in the output buffer, Qj. Q= prefix-sum(qj) = [0,(q0+q1),(q0+q1, q2), ..., 254 X j=0 qj] (3.16) Dividing these results by 32 and obtaining the remainder yields the desired word position and starting bit. Once the two coordinates are calculated, the codewords are shifted in such a way that they start in the corresponding starting bit. Afterwards, a logical OR is performed between the codewords which share the same word position, using CUDA atomic operations, which make it possible for a thread to perform an operation without interference from any other threads, to avoid having threads with the same word position accessing the output encoded buffer at the same time. 94 Chapter 3. Lossy hyperspectral compression on GPU 3.5.4 Bit packing Once the encoded blocks corresponding to a 16×16 pixel portion of the image are obtained, they have to be written to a single bit stream which represents the compressed image. As each block has been processed independently in the entropy coding kernel, the resulting encoded blocks have been written to a specific position of the global memory. To construct the final output bit stream, the encoded blocks have to be saved in sequential order and, as it happened with the codewords, in such a way that no bits are left unused between them. A similar strategy to the one used by the entropy coder is followed: the word position and starting bit is calculated for every compressed block. Afterwards, the compressed block is shifted and and atomic OR is performed. However, this time it is necessary to perform the operations on complete encoded blocks (conformed by more than one 32-bits variable). The word position and starting bit are calculated now for every encoded block. The blocks are copied from global to shared memory where they are shifted according to the starting bit in parallel. Finally the atomic operations are used to perform the logic OR and create the compressed bit stream, see Figure 3.14. The strategy to perform all these operations in parallel is explained in more detail next. Calculating the final position of a block in the output buffer Let pibe the position of the last codeword written in block i, and let libe the number of bits left unused in the 32-bits variable where the codeword in piis stored. Let NB be the total number of blocks in the image, computed as: Chapter 3. Lossy hyperspectral compression on GPU 95 Output from entropy coder (compressed 16x16 blocks) Final bit stream (content of the compressed file) OR block i block i+1 block i shifted block i+1 word 0 word 1 word piword 255 word 0 word 1 word pi+1 word 255 word 0 word 1 word pi word 0 word 1 word pi+1 word pi+1+1 lili+1 Figure 3.14: Bit packing strategy NB =Nz ×Ny ×Nx 16 ×16 (3.17) In order to calculate the number of bits that a block ihas to be shifted to the left, a prefix-sum of liis performed, obtaining for each block: Li=     0 if i= 0 i−1 P i=1 li−1if i > 0(3.18) The number of bits that a block has to be shifted to the left, sh leftiis then calculated as follows: sh lefti=   Li,if Li<32 Limod 32,if Li≥32 (3.19) where operator mod stands for modulo. The prefix-sum performed to calculate Liis computed using the same strategy explained in Section 3.5.3. 102 Chapter 3. Lossy hyperspectral compression on GPU 3.6.1 Preliminary considerations The main drawback encountered when trying to parallelize the LCE decompressor algorithm lies in the fact that the codewords have been buffered to the compressed file one after the other, in a bit-by-bit fashion and it is impossible to know when a compressed block starts. Therefore, in principle, although the blocks can be decoded independently, the compressed file must be read sequentially to obtain the codewords for every 16 ×16 block. There are several reasons why it is not feasible to find out where the compressed blocks start in the compressed file: The codewords are of variable length, therefore the compressed 16×16 blocks have different sizes in bits. A single block contains codewords obtained in different ways: –α,µzand the first sample of the first band are not encoded. –The first sample of blocks in band z > 0 is encoded using an exponential Golomb code. –The rest of the samples are encoded using a Golomb code. These facts make it impossible to split the compressed file in order to obtain the different decompressed blocks in parallel without additional information about the location of the blocks in the compressed bit stream. In order to solve this issue, it is proposed to add a header to the compressed file, which contains information about the size in bits of each compressed block. The downside of adding a header is that it increases the size of the compressed file, and therefore reduces the compression ratio. Consequently, the header must be designed in such a way that it allows achieving more parallelization, without affecting the compression ratio considerably. The amount of acceleration that can be obtained by adding the header can be calculated taking Chapter 3. Lossy hyperspectral compression on GPU 103 into account the number of blocks that will be processed in parallel and the number of iterations needed to decode each block. 3.6.2 Header design Two different options are explored before adding a header to the compressed file. For simplicity and as a worst case scenario, it is considered that the header is saved at the beginning of the compressed file, storing each value in a 32-bits word. The impact of adding a header to the compressed image is evaluated in terms of: The increment (I) in bits of the size of the compressed file, which can be calculated as: I=Size of header ×32 (3.22) The number of elements that can be processed in parallel if the header is added (PBLOCKS), which gives an idea of the amount of acceleration that will be obtained. The number of sequential iterations (SI) that are necessary to decode each parallel element. Taking into account these facts, a figure of merit (F M) is defined: FM =PBLOCKS I×SI (3.23) The figure of merit is evaluated for each of the proposed options, in order to find out which of them is able to potentially produce a higher speedup with the lowest possible impact on the compression ratio. 104 Chapter 3. Lossy hyperspectral compression on GPU Option1 In the first option explored, the header contains the size in bits of every 16 ×16 spatial block with all its bands, as shown in Figure 3.17. For this option, the size of the compressed file increases proportionally to the spatial dimension of the hyperspectral image. Hyperspectral image 16x16 block (0,0) 16x16 block (0,1) 16x16 block (1,0) Ny Nx Nz number of bits in compressed band0 number of bits in compressed band1 number of bits in compressed band2 number of bits in compressed bandNz-1 HEADER OPTION 1 COMPRESSED FILE compressed band0[block{0},...block{Nbv x Nbh-1}] compressed band1[block{0},...block{Nbv x Nbh-1}] compressed bandNz-1[block{0},...block{Nbv x Nbh-1}] Figure 3.17: Format of header and compressed file for Option1. The increment of the compressed file in this case is calculated as: I=Size of header ×32 = Ny ×Nx 16 ×16 ×32 (3.24) where Nyand Nxdenote the number of lines and columns in the hyperspectral cube respectively. Chapter 3. Lossy hyperspectral compression on GPU 105 The number of elements which can be processed in parallel corresponds to the number of 16 ×16 spatial blocks with all their bands present in the hyperspectral image. PBLOCKS =Ny ×Nx 16 ×16 (3.25) The number of sequential iterations which are necessary to decode the compressed buffer is the number of bands in the image, Nz: SI =Nz (3.26) Finally, the figure of merit (FM) is calculated, yielding the result: FMOP T 1=1 32 ×Nz (3.27) Option2 In the second option, it is proposed that the header contains the size in bits of all 16 ×16 spatial blocks in a specific band, resulting in a header with the format shown in Figure 3.18. With this header, the spatial blocks in a band can be decoded in parallel. The header contains a value for every band in the image, each saved in a 32-bits variable. Adding the described header involves making additional changes to the encoder, since the order in which the codewords are saved to the compressed file must be altered. Specifically, as it is shown in Figure 3.18, it must contain the codewords which represent a compressed band with all its spatial blocks in order. The changes must be made to the compressor’s main nested loop, which must be modified so that all horizontal and vertical 16 ×16 spatial blocks in a band are processed before the compression of the next band begins. Additional changes might be needed in order to guarantee 106 Chapter 3. Lossy hyperspectral compression on GPU that the necessary information from the previous band is available for the compression of a specific band. Hyperspectral image 16x16 block (0,0) 16x16 block (0,1) 16x16 block (1,0) Ny Nx Nz number of bits in compressed block{0} number of bits in compressed block{1} number of bits in compressed block{2} number of bits in compressed block{Nbv x Nbh-1} HEADER OPTION 2 COMPRESSED FILE compressed block{0}[band0,...bandNz-1] compressed block{1}[band0,...bandNz-1] compressed block{Nbv x Nbh-1}[band0,...bandNz-1] Figure 3.18: Format of header and compressed file for Option2. The increase of size (I) in the compressed file and the number of parallel elements (P BLOCKS) for this case are: I=Size of header ×32 = Nz ×32 (3.28) PBLOCKS =Nz (3.29) The number of sequential iterations (SI) which have to be performed by the decoder for every parallel element corresponds to the number of spatial 16 ×16 blocks, given by: Chapter 3. Lossy hyperspectral compression on GPU 107 SI =Ny ×Nx 16 ×16 (3.30) Yielding a figure of merit of: FMOP T 2=8 Ny ×Nx (3.31) Comparing the figures of merit resulting from both options, it can be observed that the best trade-off between the increase of the compression ratio and the potential speedup achieved depends basically on the image dimensions. For images which are spatially big and do not comprise a very high number of bands, which is the case of multispectral and hyperspectral images, Option1 is more convenient. However, for ultraspectral images, which can contain thousands of bands and are usually smaller in the spatial dimensions, Option2 should be considered. It is decided to adopt Option1 for the reasons stated below: Most of the sensors are multispectral or hyperspectral, therefore performance is expected to be better with Option1. Other than adding the header to the compressed file, no changes need to be made to the original LCE compressor. 3.6.3 Decoding the blocks The first stage of the LCE parallel decompressor is to read and decode the compressed file. The header containing the size in bits of each compressed 16 ×16 block with all its bands is attached to the compressed bit stream to make it possible to parallelize this stage. 108 Chapter 3. Lossy hyperspectral compression on GPU With the information in the header, the exact location of the first codeword of each 16×16 compressed block can be calculated performing a summation. A kernel is created to read and decode every compressed block in parallel in such a way that once a block is read from the compressed buffer, the kernel iterates sequentially to decode all the codewords. These iterations have to cover all bands and all samples in a 16 ×16 block. As it will be demonstrated in the experimental results (Section 3.7) this loop represents the main performance weakness of the parallel LCE decompressor. After the execution of this kernel, parameters ˆα,µzand all the decoded prediction errors are saved to the CUDA global memory so that the inverse predictor can be applied to them. 3.6.4 Inverse quantization and prediction The kernel designed to perform the inverse quantization and prediction is almost the same as the one designed for the LCE CUDA compressor, however, it performs the inverse operations in the inverse order. As for the compressor, the inverse quantization and prediction kernel is designed in such a way that a CUDA block is launched for every block in the image and every CUDA block launches 256 (16 ×16) threads. 3.7 Experimental results The CUDA implementation of the LCE compressor and decompressor are executed on an Nvidia Tesla C2075 GPU. The experiments are designed in such a way that it is possible to assess several aspects: the accuracy of the GPU implementation; the acceleration achieved when compared with the execution of the algorithm in a single threaded CPU; the relationship Chapter 3. Lossy hyperspectral compression on GPU 109 Table 3.2: Hyperspectral images under compression Sensor Area Nz Ny Nx bpppb MODIS - 17 1984 1344 12 AVIRIS Indian Pines 220 1952 608 16 AIRS Granule 1501 128 80 14 between the achieved acceleration and the hyperspectral image dimensions and the configuration parameters of the LCE algorithm. Three different hyperspectral scenes, acquired by different sensors and with different spatial and spectral sizes, are used as targets for compression. Table 5.11 summarizes the main characteristics of the images: the sensor which acquired them; the number of lines (Ny), columns (Nx) and bands (Nz); and the number of bits utilized to represent the raw samples (bpppb). All the experiments are performed in a worskstation with a 3.19 GHz Intel Xeon W5580 processor, running on a 64-bits operating system with 12 GBytes of RAM memory. The LCE algorithm implementation operates with integer variables only, except for the quantization stage, where it is necessary to utilize double precision floating point variables when the uniformthreshold quantizer (UTQ) is enabled. 3.7.1 Validation Before any experiment is performed, it is mandatory to verify that the GPU implementations of both the compressor and decompressor produce the same results as their CPU version counterpart. For this purpose, the images are compressed and decompressed with the GPU and the CPU implementations. The resulting compressed and reconstructed files are then compared bit by bit, demonstrating that the results of both implementations are identical. 110 Chapter 3. Lossy hyperspectral compression on GPU 3.7.2 Impact of adding a header Once it is verified that the GPU versions of the LCE algorithm produce the correct results, we assess the impact in the compression ratio of adding a header to the compressed file, in order to demonstrate the validity of the decisions adopted in Section 3.6.2. The compression strength of the LCE algorithm can be configured by the user by setting a parameter delta. Low values of delta yield low compression while a high delta value increases the compression ratio at the cost of producing more losses of information. The impact of adding a header is more significant if the compression ratio is high, i.e. the compressed file is smaller, hence, each hyperspectral image under test is compressed with two different values of delta, namely 1 and 60 in order to perform a fair evaluation. The compression ratio (CR) in bits per pixel per band (bpppb) is calculated as: CR (bpppb) = Size of compressed file (bits) Nz ×Ny ×Nx (3.32) Moreover, the percentage of increment in the compression ratio (CR increment) is calculated as: CR increment(%) = Size of header (bits) Size of compressed file (bits)×100 (3.33) It is observed in Table 3.3 that, although adding a header to the compressed file increases the compression ratio, the impact produced is not considerable, with an average of 0.06% and a maximum of 0.20% for all the images under test. Chapter 3. Lossy hyperspectral compression on GPU 111 Table 3.3: Impact of adding a header to the compressed file Sensor Size (Nz ×Ny ×Nx) delta CR without header CR with header Increment (%) MODIS 17 ×1984 ×1344 1 7.3352 7.3426 0.10 60 3.9538 3.9612 0.19 AVIRIS 220 ×1952 ×608 1 5.8970 5.8976 0.01 60 2.3604 2.3610 0.02 AIRS 1501 ×125 ×80 1 4.3474 4.3475 0.00 60 0.8704 0.8704 0.01 3.7.3 Profiling The GPU implementation of the LCE compressor and decompressor are profiled with the software tools supplied by Nvidia. For this purpose a 512× 512 subimage of the AVIRIS scene, which comprises 220 bands, is employed. This hyperspectral cube is compressed with the GPU implementations and profiled with the tools supplied by Nvidia, in order to detect which of the kernels is the most time-consuming and where are the main bottlenecks of the implementation. GPU compressor profiling Figure 3.19 shows the profiling of the LCE compressor. The total time to perform the compression is 394 ms, which is almost a quarter of the time achieved with the GPU implementation of JPEG2000 presented in [19]. The predictor takes a 23.29%, what shows that, although the predictor has to loop to cover all bands in the hyperspectral cube, this fact is not a bottleneck for global performance. The most time-consuming operations for both GPUs are the entropy coding (44.24%), and the memory transactions between the CPU and the GPU. The strategy designed for the bit packing, which is used to create the final bit stream shows very good results, taking only 4.73% of the total compression time. 214 Appendix A. Sinopsis en espa˜nol hardware disponible a bordo de los sat´elites: alta eficiencia en la compresi´on y baja complejidad. Existen asimismo otros algoritmos, espec´ıficamente dise˜nados para la compresi´on de im´agenes hiperespectrales a bordo de sat´elites, en la literatura cient´ıfica disponible actualmente. Entre ellos, cabe destacar por su relevancia en el trabajo presentado en esta Tesis Doctoral, el algoritmo para compresi´on con p´erdidas conocido como Lossy Compression for Exomars (LCE). Este algoritmo est´a basado en t´ecnicas de predici´on, e introduce p´erdidas de informaci´on combinando la t´ecnica de cuantificaci´on y una t´ecnica que ajusta el ratio de compresi´on. Este ajuste permite que la distorsi´on en la imagen reconstruida no exceda un determinado umbral seleccionado por el usuario. No solamente es importante desarrollar nuevos algoritmos para la compresi´on de im´agenes hiperespectrales a bordo de sat´elites. Tambi´en es necesario considerar el soporte f´ısico (hardware) en que se van a ejecutar, ya que debe cumplir con ciertos requisitos, entre los que destacan la tolerancia a la radiaci´on solar y un bajo consumo de potencia. Generalmente, los algoritmos se implementan en procesadores de a bordo (on-board processor), procesadores digitales de se˜nal (digital signal processor, DSP), circuitos integrados de aplicaci´on espec´ıfica (application-specific integrated circuit, ASIC), o field programable gate array (FPGA). Una FPGA es un dispositivo que contiene bloques de l´ogica programable cuya interconexi´on y funcionalidad puede ser configurada. Recientemente, las FPGAs han alcanzado gran popularidad para la implementaci´on de aplicaciones en el sector aeroespacial, principalmente debido a su alto rendimiento, bajo consumo y la posibilidad de ser reprogramadas. Adem´as, existen versiones en el mercado tolerantes a la radiaci´on y, por lo tanto, aptas para operar en un sat´elite. Aparte de los soportes f´ısicos previamente enumerados, en el ´ambito de la computaci´on de altas prestaciones se ha popularizado recientemente la Appendix A. Sinopsis en espa˜nol 215 programaci´on de prop´osito general en tarjetas gr´aficas (graphics processing unit (GPU)), conocidas como GPU por sus siglas en ingl´es. Inicialmente, las GPUs se hicieron populares en la industria de los videojuegos por su alta capacidad para el procesamiento masivo de datos en paralelo. En los ´ultimos a˜nos se ha ido extendiendo su uso hacia la computaci´on de prop´osito general, gracias en parte a la aparici´on de entornos y directivas para la programaci´on de estos dispositivos, como es el caso de Computer Unified Device Architecture (CUDA), desarrollado para la programaci´on de tarjetas gr´aficas de NVidia. Una GPU consiste en un conjunto de multiprocesadores que trabajan en paralelo en modo SIMD (del ingl´es Single Instruction, Multiple Data, en espa˜nol: “una instrucci´on, m´ultiples datos”). A diferencia de los procesadores de pr´oposito general o CPUs, las GPUs dedican la mayor parte de sus componentes al c´omputo en vez de al control y la memoria, por lo que consiguen acelerar sustancialmente los algoritmos caracterizados por su alto paralelismo a nivel de datos. Tal es el caso de los algoritmos de compresi´on de im´agenes hiperespectrales, que bien podr´ıan utilizar esta tecnolog´ıa para reducir sus tiempos de c´omputo. En este sentido, uno de los objetivos de este trabajo de Tesis ser´a estudiar c´omo de convenientes son las GPUs para acelerar los algoritmos de compresi´on de im´agenes hiperespectrales. A pesar de que actualmente las GPUs no est´an cualificadas para trabajar en el espacio, debido a su elevado consumo de potencia y por no ser tolerantes a la radiaci´on solar, tienen un elevado potencial. Por lo tanto, no se descarta que en un futuro se lleven a cabo los desarrollos tecnol´ogicos necesarios para habilitarlas para ser utilizadas para la compresi´on de datos a bordo de sat´elites en misiones espaciales. A.2 Objetivos y metodolg´ıa de trabajo Se espera que los sensores hiperespectrales en el futuro capturen cada vez m´as volumen de datos, mientras que los anchos de banda de transmisi´on y la 216 Appendix A. Sinopsis en espa˜nol cantidad de almacenamiento disponible se mantendr´an relativamente estables. Como resultado, la compresi´on de datos a bordo se har´a indispensable. En este sentido, una compresi´on eficiente depender´a, tanto del algoritmo utilizado, como de la tecnolog´ıa f´ısica en que se ejecute dicho algoritmo. La meta principal de este trabajo de Tesis es aportar nuevas soluciones para la implementaci´on f´ısica y ejecuci´on de algoritmos de im´agenes hiperespectrales a bordo de sat´elites. Este trabajo de investigaci´on ayudar´a a mejorar los actuales resultados de las implemetaciones hardware existentes dentro del estado del arte, estableciendo adem´as directrices para los futuros trabajos de investigaci´on en este ´ambito. A continuaci´on se enumeran los objetivos de esta Tesis: Caracterizar los algoritmos de compresi´on de im´agenes hiperespectrales. Proponer soluciones para acelerar dichos algoritmos en distintas tecnolog´ıas hardware, en concreto GPUs y FPGAs. Comparar las soluciones en ambas tecnolog´ıas en cuanto a su rendimiento, coste de la soluci´on, flexibilidad de la implementaci´on y consumo de potencia. Validar las soluciones propuestas, mostrando el correcto funcionamiento de las implementaciones resultantes. Para la consecuci´on de estos objetivos se realiza la implementaci´on de dos algoritmos, dise˜nados espec´ıficamente para la compresi´on de im´agenes hiperespectrales a bordo de sat´elites, en dos tecnolog´ıas hardware: FPGAs y GPUs. Ambas tecnolog´ıas son capaces de acelerar algoritmos que permiten el procesado de datos en paralelo, pero presentan m´ultiples diferencias en cuanto a sus metodolog´ıas de dise˜no, consumo de potencia y tolerancia a la radiaci´on solar. En concreto, se dise˜na la implementaci´on en GPU y FPGA Appendix A. Sinopsis en espa˜nol 217 de un algoritmo de compresi´on sin p´erdidas de im´agenes hiperespectrales, lo que permite hacer una comparativa en t´erminos de rendimiento entre ambas alternativas tecnol´ogicas. Finalmente, se realiza la implementaci´on del algoritmo est´andar, el CCSDS para compresi´on sin p´erdidas de im´agenes hiperespectrales, en una FPGA cualificada para el espacio. Este trabajo sirve para evaluar la viabilidad de la utilizaci´on del est´andar CCSDS 123 en misiones espaciales actuales, y es a su vez, sirve como ejemplo de caso de uso. Se realizar´a una comparativa entre los dos algoritmos estudiados y las tecnolog´ıas hardware consideradas, as´ı como con otras soluciones del estado del arte. A.3 Compresi´on con p´erdidas en GPU y FPGA Con el objetivo de estudiar las posibles implementaciones hardware para la compresi´on de im´agenes hiperespectrales, a bordo de sat´elites, se proponen en este trabajo diferentes soluciones. En concreto, se realiza la implementaci´on del algoritmo LCE en una GPU y una FPGA. El algoritmo LCE es particularmente relevante, ya que fue dise˜nado espec´ıficamente para ser ejecutado a bordo de sat´elites. Por lo tanto, cumple con determinados requisitos que son cruciales para poder realizar la compresi´on a bordo de un sat´elite, como son: baja complejidad, alta eficiencia, tolerancia a fallos y sencillez de cara a la implementaci´on en dispositivos hardware. El algoritmo comprime la imagen por bloques independientes de 16 ×16 p´ıxeles con todas sus bandas. Consta de una primera etapa de predicci´on en la que se obtienen los errores de predicci´on, los cuales se codifican mediante un codificador entr´opico. Las palabras de c´odigo obtenidas se empaquetan en un fichero que representa la imagen comprimida. La t´ecnica de compresi´on por bloques hace que el LCE sea tolerante a fallos, ya que un error en la transmisi´on de un bloque no impide la decodificaci´on del resto de 218 Appendix A. Sinopsis en espa˜nol bloques, una vez la imagen sea recibida. La divisi´on de la imagen en bloques se hace como se muestra en la Figura A.1. 16x16 block (0,0) 16x16 block (0,1) 16x16 block (1,0) sz,y,x Ny Nx Nz z y x Nz z m n 16 16 sz,m,n Hyperspectral image Block Figure A.1: Divisi´on de la imagen hiperespectral en bloques independientes Al inicio de este trabajo, se contaba con una implementaci´on software del algoritmo LCE descrito en lenguaje de programaci´on C, que es posible ejecutar en cualquier PC o estaci´on de trabajo. Este software sirve como base para desarrollar tanto la implementaci´o en GPU como la implementaci´on en FPGA, y adem´as se utiliza como referencia para validar los resultados obtenidos con ambas tecnolog´ıas. Appendix A. Sinopsis en espa˜nol 219 A.3.1 Implementaci´on del algoritmo LCE en una GPU Se realiza una evaluaci´on de la viabilidad de las GPUs para ejecutar la compresi´on de im´agenes hiperespectrales, a trav´es de una paralelizaci´on y ejecuci´on del algoritmo LCE en una GPU de Nvidia utilizando CUDA . Durante la paralelizaci´on, se adapta el algoritmo a las abstracciones de CUDA, dise˜nando un kernel paralelo para la ejecuci´on de cada una de las etapas del algoritmo. Esto implica aplicar diferentes estrategias para identificar qu´e datos pueden procesarse independientemente, y qu´e dependencias de datos es posible eliminar. Adem´as, se utilizan estrategias para la reduci´on de iteraciones en los bucles en que las dependencias de datos son inevitables. Todo esto se realiza teniendo en cuenta la arquitectura hardware de la GPU seleccionada, de modo que pueda maximizarse la utilizaci´on de los multiprocesadores presentes en ella. La implementaci´on dise˜nada muestra una alta aceleraci´on, hasta 15.41 veces m´as r´apida que la ejecuci´on del mismo algoritmo en una CPU. El rendimiento, definido como el n´umero de muestras procesadas por unidad de tiempo, es igualmente muy alto. El algoritmo LCE es capaz de comprimir entre 100 y 140 muestras cada segundo cuando se ejecuta en una GPU. Las dificultades y cuellos de botella m´as importantes encontrados a la hora de realizar la paralelizaci´on del LCE se han identificado como parte de este trabajo. Destaca principalmente la gran cantidad de tiempo empleado en el env´ıo de la imagen hiperespectral sin comprimir a la GPU para que pueda ser procesada. Se proponen l´ıneas de trabajo futuras para superar este problema y optimizar el rendimiento de la implementaci´on. Teniendo en cuenta las similitudes entre el compresor LCE y el correspondiente decompresor, se considera conveniente relizar tambi´en una implementaci´on en GPU de dicho decompresor. La paralelizaci´on del decompresor 220 Appendix A. Sinopsis en espa˜nol supone dificultades adicionales, que hacen que el rendimiento obtenido por la GPU sea m´as bajo que el que se obtuvo para el compresor. En concreto, el principal problema con que se encuentra el decompresor es el formato de la imagen comprimida, que impide la identificaci´on de los elementos independientes que es necesario extraer para poder realizar las tareas de decompresi´on en paralelo. En este trabajo de Tesis se propone una soluci´on parcial a este problema, que consiste en a˜nadir una cabecera al archivo comprimido con informaci´on relativa a la ubicaci´on de los elementos independientes en el archivo comprimido. La inclusi´on de esta cabecera permite procesar m´as datos en paralelo, pero a su vez aumenta el tama˜no del archivo comprimido, reduciendo por tanto el ratio de compresi´on. Se estudian distintos formatos de cabecera, eligiendo aquel que permite obtener la mejor relaci´on entre la cantidad de elementos que se descomprimen en paralelo y el tama˜no de la cabecera. Como ya se ha comentado anteriormente, hoy en d´ıa no es posible utilizar las GPU para la compresi´on a bordo de sat´elites, ya que no son tolerantes a la radiaci´on y presentan consumos de potencia demasiado elevados. Sin embargo, este trabajo de investigaci´on es ´util pra demonstrar que es posible obtener implementaciones eficientes y flexibles cuando se utilizan las GPUs para la compresi´on de datos de im´agenes hiperspectrales, y contribuye a justificar y motivar la investigaci´on en este campo, de manera que se hagan esfuerzos para desarrollar GPUs cualificadas para el espacio en un futuro. A.3.2 Implementaci´on del algoritmo LCE en una FPGA Se realiza adem´as una implementaci´on del algoritmo LCE para compresi´on de im´agenes hiperespectrales con p´erdidas a bordo de sat´elites sobre una FPGA. Para realizar dicha implementaci´on se utiliza la herramienta de Appendix A. Sinopsis en espa˜nol 221 s´ıntesis autom´atica CatapultC, que genera los ficheros de descripci´on hardware a partir de c´odigo escrito en lenguaje de programaci´on C. Esto hace posible obtener implementaciones en un tiempo reducido. La metodolog´ıa de trabajo en CatapultC consta de varias etapas. La primera de ellas es la optimizaci´on y verificaci´on del c´odigo en lenguaje C. Seguidamente se establecen las restricciones globales del hardware, indicando la frecuencia de reloj, el comportamiento de la se˜nal de reset y la tecnolog´ıa en que se va a realizar la implementaci´on final. Se debe indicar en este paso qu´e funci´on dentro del c´odigo C representa el bloque m´as alto de la jerarqu´ıa. A continuaci´on se establecen restricciones concretas de la arquitectura, realiz´andose, entre otras tareas, la optimizaci´on de bucles, que pueden ser desenrrollados (unrolling) o segmentados (pipeline). Finalmente, CatapultC realiza la planificaci´on temporal del dise˜no y genera la descripci´on en VHDL del c´odigo. ´ Esta es v´alidada antes de ser ser introducida en herramientras de s´ıntesis, como Mentor Graphis Precision o Synplify. Estas herramientas generan los archivos necesarios para programar la FPGA con el dise˜no final. Las modificaciones realizadas al c´odigo C original del compresor LCE fueron tanto sint´acticas como funcionales. Las modificaciones sint´acticas se centran en eliminar aquellos elementos del lenguaje C no reconocibles por la herramienta CatapultC, ya que no son realizables en hardware; como puede ser el alojamiento de memoria de forma din´amica. Tambi´en se realizan modificaciones funcionales, cambiando el algoritmo para evitar realizar computaciones que en hardware requieren un alto consumo de recursos, como puede ser una divisi´on de dos variables enteras. Una vez realizadas las modificaciones en el c´odigo C, se procede a optimizar la configuraci´on en la herramienta CatapultC para obtener el resultado deseado en t´erminos de ´area y latencia. El m´odulo implementado realiza la compresi´on de un bloque de 16 ×16 muestras. Sus entradas y salidas se representan en la Figura A.2. Adem´as 222 Appendix A. Sinopsis en espa˜nol Table A.1: Ocupaci´on en la FPGA del LCE Virtex 5VFX130 Modular approach Non-modular approach Resources Available Used %Used % DSP48Es 320 17 5 25 8 Number of RAMB18X2s 298 4 1 4 1 Number of slices 20480 2951 14 1935 10 Number of Slice Registers 81920 4208 5 5995 7 Number of Slice LUTS 81920 7836 9 7738 10 de la implementaci´on de este m´odulo con CatapultC, se realiza una implementaci´on modular del mismo algoritmo LCE, generando la descripci´on VHDL de los distintos bloques funcionales con CatapultC y dise˜nando el control y flujo de datos entre los m´odulos manualmente con el objetivo de mejorar el rendimiento del dise˜no. LCE compressor ym filecountfilecount pp mm block_out RAM ref_block RAM cur_block RAM Figure A.2: M´odulo de compresi´on y sus interfaces de entrada/salida La ocupaci´on de recursos, frecuencia m´axima y rendimiento de la implementaci´on del LCE en la FPGA, utilizado las mencionadas estrategias, se muestran en las Tablas A.1 y A.2 Appendix A. Sinopsis en espa˜nol 223 Table A.2: Muestras comprimidas por segundo de la implementaci´on sobre FPGA del LCE Virtex 5VFX130 Dise˜no modular Dise˜no no modular Max.Frecuencia (MHz) 86 80 Mega-Muestras/seg 27.7 16.7 Los resultados demuestran la validez del algoritmo LCE para compresi´on a bordo de sat´elites, y prueban que es posible implementarlo con una baja utilizaci´on de recursos y baja potencia en una FPGA. Asimismo, puede observarse que, la implementaci´on modular ha conseguido mejorar el rendimiento del dise˜no, manteniendo resultados similares de ocupaci´on de la FPGA. Se realiza finalmente una comparativa de las implementaciones realizadas del compresor LCE en t´erminos de rendmiento (n´umero de muestras procesadas por unidad de tiempo). Los resultados correspondientes a la compresi´on de una imagen del sensor AVIRIS de 220 bandas mediante el algoritmo LCE sobre CPU, GPU y FPGA se muestran en la Figura A.3. En la citada figura se muestran el n´umero de muestras procesadas por segundo en funci´on del tama˜no espacial de la im´agen (n´umero de l´ıneas por n´umero de columnas). Aunque la GPU es capaz de procesar mayor n´umero de muestras por unidad de tiempo, actualmente no es posible utilizarlas en el espacio y, por lo tanto, de momento no pueden considerarse una alternativa viable a las FPGA. Adem´as, es posible acelerar a´un m´as la implementaci´on presentada en la FPGA, realizando la implementaci´on de m´as de un m´odulo compresor en el mismo dispositivo, y haciendo que trabajen en paralelo. De este modo, podr´ıa conseguirse un rendimiento comparable al obtenido en la GPU, mientras que el consumo de potencia ser´ıa notablemente inferior en la FPGA. 230 Appendix B. Publications B.1 Journals [1] Santos, L., Lopez, S., Callico, G. M., L´opez, J. F., and Sarmiento, R. (2012). Performance Evaluation of the H.264/AVC Video Coding Standard for Lossy Hyperspectral Image Compression. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, 5(2), 451–461. [2] Santos, L., Magli, E., Vitulli, R., L´opez, J. F., and Sarmiento, R. (2013). Highly-Parallel GPU Architecture for Lossy Hyperspectral Image Compression. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of,6(2), 670–681. [3] Santos, L., Magli, E., Vitulli, R., N´u˜nez, A., L´opez, J. F., and Sarmiento, R. (2013). Lossy hyperspectral image compression on a graphics processing unit: parallelization strategy and performance evaluation. Journal of Applied Remote Sensing, 7(1), 074599-074599. B.2 International Conferences [1] Santos, L., L´opez, S., Callic´o, G. M., L´opez, J. F., and Sarmiento, R. (2011). Hyperspectral image compression with a H.264/AVC baseline encoder. In Conference on Design of Circuits and Integrated Systems (DCIS 2011). [2] Santos, L., L´opez, S., Callic´o, G. M., L´opez, J. F., and Sarmiento, R. (2011). Lossy hyperspectral image compression with state-of-the-art video encoder. In Proceedings of SPIE (Vol. 8183, 81830). Appendix B. Publications 231 [3] Santos, L., Vitulli, R., L´opez, J. F., and Sarmiento, R. (2012). GPU implementarion of a lossy compression algorithm for hyperspectral images. In IEEE Workshop on Hyperspectral Image and Signal Processing – Evolution in Remote Sensing (WHISPERS) 2012. [4] Santos, L., Vitulli, R., L´opez, J. F., and Sarmiento, R. (2012). CUDA based GPU implemenmtation of a parallelized algorithm for lossy hyperspectral image compression. In On Board Payload Data Compression Workshop – OBPDC 2012. [5] Santos, L., Vitulli, R., L´opez, S., Marrero, G. M., L´opez, J. F., and Sarmiento, R. (2012). Accelerating lossy hyperspectral image compression on a GPU. In Conference on Design of Circuits and Integrated Systems (DCIS 2012). [6] Santos, L., Lopez, J. F., Sarmiento, R., and Vitulli, R. (2013). FPGA implementation of a lossy compression algorithm for hyperspectral images with a high-level synthesis tool. In Adaptive Hardware and Systems (AHS), 2013 NASA/ESA Conference on (pp. 107–114). [7] Garc´ıa, A., Santos, L., L´opez, S., Callic´o, G.M., L´opez, J.F., Sarmiento, R. (2013) High level modular implementation of a lossy hyperspectral image compression algorithm on a FPGA. In IEEE Workshop on Hyperspectral Image and Signal Processing Evolution in Remote Sensing (WHISPERS), 2013 [8] Garc´ıa, A., Santos, L., L´opez,, S., Callic´o, G. M., Lopez, J. F., and Sarmiento, R. (2014). Efficient lossy compression implementations of hyperspectral images: tools, hardware platforms, and comparisons. In SPIE Sensing Technology+ Applications (pp. 912408-912408) References [1] National Aeronautics and Space Administration (NASA). Landsat science, 2014. URL http://landsat.gsfc.nasa.gov/. [2] Centre National d’Etudes Spatiales (CNES). Spot earth observation mission. eyes in space., 2013. URL http://smsc.cnes.fr/SPOT/. [3] Indian Space Research Organisation (ISRO). Indian remote sensing (irs) satellite system, 2008. URL http://www.isro.org/scripts/ currentprogrammein.aspx#IRS. [4] National Aeronautics and Space Administration (NASA). Airborne visible/infrared imaging spectometer (aviris), 2014. URL http:// aviris.jpl.nasa.gov/. [5] Alexander FH Goetz. Three decades of hyperspectral remote sensing of the earth: A personal view. Remote Sensing of Environment, 113: S5–S16, 2009. [6] Peg Shippert. Introduction to hyperspectral image analysis. Online Journal of Space Communication, 3, 2003. URL http:// spacejournal.ohio.edu/pdf/shippert.pdf. [7] William F Belokon, Spectral Imagery Training Center, and Logicon Geodynamics. Multispectral imagery reference guide. Logicon Geodynamics, 1997. 233 234 References [8] Centre National d’Etudes Spatiales (CNES). URL http://smsc. cnes.fr/PLEIADES/. [9] European Space Agency (ESA). Proba missions, 2014. URL http: //www.esa.int/Our_Activities/Technology/Proba_Missions. [10] J Bioucas-Dias, Antonio Plaza, Gustavo Camps-Valls, Paul Scheunders, N Nasrabadi, and JOCELYN Chanussot. Hyperspectral remote sensing data analysis and future challenges. Geoscience and Remote Sensing Magazine, IEEE, 1(2):6–36, 2013. [11] Antonio Plaza, Jon Atli Benediktsson, Joseph W Boardman, Jason Brazile, Lorenzo Bruzzone, Gustavo Camps-Valls, Jocelyn Chanussot, Mathieu Fauvel, Paolo Gamba, Anthony Gualtieri, et al. Recent advances in techniques for hyperspectral image processing. Remote Sensing of Environment, 113:S110–S122, 2009. [12] J.M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Qian Du, P. Gader, and J. Chanussot. Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, 5(2):354–379, April 2012. ISSN 1939-1404. doi: 10.1109/JSTARS.2012.2194696. [13] N. Keshava and J.F. Mustard. Spectral unmixing. Signal Processing Magazine, IEEE, 19(1):44–57, Jan 2002. ISSN 1053-5888. doi: 10. 1109/79.974727. [14] David A Landgrebe. Signal Theory Methods in Multispectral Remote Sensing. Wiley, Newark, NJ, 2005. [15] G. Shaw and D. Manolakis. Signal processing for hyperspectral image exploitation. Signal Processing Magazine, IEEE, 19(1):12–16, Jan 2002. ISSN 1053-5888. doi: 10.1109/79.974715. References 235 [16] Stefania Matteoli, Marco Diani, and Giovanni Corsini. A tutorial overview of anomaly detection in hyperspectral images. Aerospace and Electronic Systems Magazine, IEEE, 25(7):5–28, 2010. [17] David WJ Stein, Scott G Beaven, Lawrence E Hoff, Edwin M Winter, Alan P Schaum, and Alan D Stocker. Anomaly detection from hyperspectral imagery. Signal Processing Magazine, IEEE, 19(1):58–69, 2002. [18] Dimitris Manolakis and Gary Shaw. Detection algorithms for hyperspectral imaging applications. Signal Processing Magazine, IEEE, 19 (1):29–43, 2002. [19] A. Plaza M. Ciznicki, K. Kurowski. Gpu implementation of jpeg2000 for hyperspectral image compression. In SPIE - The International Society for Optical Engineering, volume 8183, 2011. [20] L. Santos, S. Lopez, G.M. Callico, J.F. Lopez, and R. Sarmiento. Performance evaluation of the h.264/avc video coding standard for lossy hyperspectral image compression. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, 5(2):451 –461, april 2012. ISSN 1939-1404. doi: 10.1109/JSTARS.2011.2173906. [21] D.A. Huffman. A method for the construction of minimum-redundancy codes. Proceedings of the IRE, 40(9):1098–1101, Sept 1952. ISSN 00968390. doi: 10.1109/JRPROC.1952.273898. [22] S. Golomb. Run-length encodings (corresp.). Information Theory, IEEE Transactions on, 12(3):399–401, Jul 1966. ISSN 0018-9448. doi: 10.1109/TIT.1966.1053907. [23] R. Rice and J. Plaunt. Adaptive variable-length coding for efficient compression of spacecraft television data. Communication Technology, IEEE Transactions on, 19(6):889–897, December 1971. ISSN 00189332. doi: 10.1109/TCOM.1971.1090789. 236 References [24] Lossless data compression recommended standard CCSDS 121.0-B-2. The Consultative Committee for Space Data Systems, 2012. [25] Lossless multispectral and hyperspectral image compression recommmended standard CCSDS 123.0-B-1. The Consultative Committee for Space Data Systems, 2011. [26] C. E. Shannon. A mathematical theory of communication. SIGMOBILE Mob. Comput. Commun. Rev., 5(1):3–55, January 2001. ISSN 1559-1662. doi: 10.1145/584091.584093. URL http://doi.acm.org/ 10.1145/584091.584093. [27] W Schober, F Lansing, K Wilson, and E Webb. High data rate instrument study. JPL publication, pages 99–4, 1999. [28] R. Trautner. Esa’s roadmap for next generation payload data processors. In Proceedings of DASIA 2011 Conference, 2011. [29] European Space Agency (ESA). Euclid mission, 2014. URL http: //sci.esa.int/euclid/. [30] A.M. Di Giorgio, S.J. Liu, G. Giusi, and G. Palamara. Euclid visible imager on-board lossless data compression: Performance assessment trade-off activities. In Proceedings of 2012 ESA workshop on Onboard Payload Data Compression (OBPDC), 2012. [31] C.M. Hartzell, L.C. Graham, T.S. Tao, H.R. Goldberg, J. CarpenaNunez, D.M. Racek, C.E. Taylor, and C.D. Norton. Data system design for a hyperspectral imaging mission concept. In Aerospace conference, 2009 IEEE, pages 1–21, March 2009. doi: 10.1109/AERO.2009. 4839507. [32] N. Aranki, A. Bakhshi, D. Keymeulen, and M. Klimesh. Fast and adaptive lossless on-board hyperspectral data compression system for References 237 space applications. In Aerospace conference, 2009 IEEE, pages 1–8, 2009. doi: 10.1109/AERO.2009.4839534. [33] A. Abrardo, M. Barni, and E. Magli. Low-complexity predictive lossy compression of hyperspectral and ultraspectral images. In Acoustics, Speech and Signal Processing (ICASSP), 2011 IEEE International Conference on, pages 797 –800, may 2011. doi: 10.1109/ICASSP.2011. 5946524. [34] Image data compression recommended standard CCSDS 122.0-B-1. The Consultative Committee for Space Data Systems, 2005. [35] M.J. Weinberger, G. Seroussi, and G. Sapiro. The loco-i lossless image compression algorithm: principles and standardization into jpeg-ls. Image Processing, IEEE Transactions on, 9(8):1309–1324, 2000. ISSN 1057-7149. doi: 10.1109/83.855427. [36] Xiaolin Wu and N. Memon. Context-based, adaptive, lossless image coding. Communications, IEEE Transactions on, 45(4):437–444, 1997. ISSN 0090-6778. doi: 10.1109/26.585919. [37] S. Hunt and L.S. Rodriguez. Fast piecewise linear predictors for lossless compression of hyperspectral imagery. In Geoscience and Remote Sensing Symposium, 2004. IGARSS ’04. Proceedings. 2004 IEEE International, volume 1, pages –312, 2004. doi: 10.1109/IGARSS.2004. 1369023. [38] Xiaolin Wu and N. Memon. Context-based lossless interband compression-extending calic. Image Processing, IEEE Transactions on, 9(6):994–1001, 2000. ISSN 1057-7149. doi: 10.1109/83.846242. [39] B. Penna, T. Tillo, E. Magli, and G. Olmo. Transform coding techniques for lossy hyperspectral data compression. Geoscience and Remote Sensing, IEEE Transactions on, 45(5):1408 –1421, may 2007. ISSN 0196-2892. doi: 10.1109/TGRS.2007.894565. 238 References [40] Bruno Aiazzi, Luciano Alparone, and Stefano Baronti. Quality issues for compression of hyperspectral imagery through spectrally adaptive dpcm. In Bormin Huang, editor, Satellite Data Compression, pages 115–147. Springer New York, 2011. ISBN 978-1-4614-11826. doi: 10.1007/978-1-4614-1183-3 6. URL http://dx.doi.org/10. 1007/978-1-4614-1183-3_6. [41] Chulhee Lee, Sangwook Lee, and Jonghwa Lee. Effects of lossy compression on hyperspectral classification. In Bormin Huang, editor, Satellite Data Compression, pages 269–285. Springer New York, 2011. ISBN 978-1-4614-1182-6. doi: 10.1007/978-1-4614-1183-3 13. URL http://dx.doi.org/10.1007/978-1-4614-1183-3_13. [42] F. Garcia-Vilchez, J. Mu˜noz-Mari, M. Zortea, I. Blanes, V. GonzalezRuiz, G. Camps-Valls, A. Plaza, and J. Serra-Sagrista. On the impact of lossy compression on hyperspectral image classification and unmixing. Geoscience and Remote Sensing Letters, IEEE, 8(2):253– 257, 2011. ISSN 1545-598X. doi: 10.1109/LGRS.2010.2062484. [43] Q. Du, N. Ly, and J.E. Fowler. An operational approach to pca+jpeg2000 compression of hyperspectral imagery. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, PP(99):1–9, 2013. ISSN 1939-1404. doi: 10.1109/JSTARS.2013. 2274527. [44] Jarno Mielikainen. Lookup-table based hyperspectral data compression. In Bormin Huang, editor, Satellite Data Compression, pages 169–184. Springer New York, 2011. ISBN 978-1-4614-1182-6. doi: 10.1007/978-1-4614-1183-3 8. URL http://dx.doi.org/10.1007/ 978-1-4614-1183-3_8. [45] D.S. Taubman and M.W. Marcellin. JPEG2000: Image Compression Fundamentals, Standards, and Practice. Kluwer, 2001. References 239 [46] L. Chang, Ching-Min Cheng, and Ting-Chung Chen. An efficient adaptive klt for multispectral image compression. In Image Analysis and Interpretation, 2000. Proceedings. 4th IEEE Southwest Symposium, pages 252–255, 2000. doi: 10.1109/IAI.2000.839610. [47] Pengwei Hao and Q. Shi. Reversible integer klt for progressive-tolossless compression of multiple component images. In Image Processing, 2003. ICIP 2003. Proceedings. 2003 International Conference on, volume 1, pages I–633–6 vol.1, 2003. doi: 10.1109/ICIP.2003.1247041. [48] B. Penna, T. Tillo, E. Magli, and G. Olmo. Progressive 3-d coding of hyperspectral images based on jpeg 2000. Geoscience and Remote Sensing Letters, IEEE, 3(1):125–129, 2006. ISSN 1545-598X. doi: 10.1109/LGRS.2005.859942. [49] N.R.M. Noor and T. Vladimirova. Parallelised fault-tolerant integer klt implementation for lossless hyperspectral image compression on board satellites. In Adaptive Hardware and Systems (AHS), 2013 NASA/ESA Conference on, pages 115–122, 2013. doi: 10.1109/AHS.2013.6604234. [50] Qian Du and J.E. Fowler. Hyperspectral image compression using jpeg2000 and principal component analysis. Geoscience and Remote Sensing Letters, IEEE, 4(2):201 –205, april 2007. ISSN 1545-598X. doi: 10.1109/LGRS.2006.888109. [51] Wei Zhu, Qian Du, and JamesE. Fowler. Hyperspectral image compression using segmented principal component analysis. In Bormin Huang, editor, Satellite Data Compression, pages 233–251. Springer New York, 2011. ISBN 978-1-4614-1182-6. doi: 10.1007/978-1-4614-1183-3 11. URL http://dx.doi.org/10.1007/978-1-4614-1183-3_11. [52] A. Karami, M. Yazdi, and G. Mercier. Compression of hyperspectral images using discerete wavelet transform and tucker decomposition. Selected Topics in Applied Earth Observations and Remote Sensing, 246 References [96] Nvidia Corporation. Nvidia’s next generation cuda compute architecture fermi, 2009. URL http://www.nvidia.com/content/ PDF/fermi_white_papers/NVIDIA_Fermi_Compute_Architecture_ Whitepaper.pdf. [97] M. Harris. Optimizing cuda. Supercomputing conference, Reno, NV, 2007. [98] M. Harris and M. Garland. Optimizing parallel prefix operations for the fermi architecture, 2011. [99] A. Balevic. Parallel variable-length encoding on gpgpus. In Proceedings of Euro-Par’09 Int. Conf. on Parallel Processing, pages 26–35, 2009. [100] Mark Harris, Shubhabrata Sengupta, and John D Owens. Parallel prefix sum (scan) with cuda. GPU gems, 3(39):851–876, 2007. [101] Catapult R C Synthesis User’s and Reference Manual - University Version - Release 2010a. Mentor Graphics Corporation, 2010. [102] Sangchul Kim, Hyunjin Kim, Taeil Chung, and Jin-Gyeong Kim. Design of h.264 video encoder with c to rtl design tool. In SoC Design Conference (ISOCC), 2012 International, pages 171–174, 2012. doi: 10.1109/ISOCC.2012.6407067. [103] T. Damak, I. Werda, N. Masmoudi, and S. Bilavarn. Fast prototyping h.264 deblocking filter using esl tools. In Systems, Signals and Devices (SSD), 2011 8th International Multi-Conference on, pages 1–4, 2011. doi: 10.1109/SSD.2011.5767375. [104] J.E. Sanchez, Aug´e E., Kiely A., Blanes I., and J. Serra-Sagrist`a. Performance impact of parameter tuning on the ccsds-123 lossless multiand hyperspectral image compression standard. In Proceedings of 2012 ESA workshop on Onboard Payload Data Compression (OBPDC), 2012. References 247 [105] Group on Interactive Coding of Images Universitat Autonoma de Barcelona. Emporda software (a ccsds-123 implementation), 2013. URL http://gici.uab.es/GiciWebPage/emporda.php. [106] Eumetsat. Eumetsat polar system - second generation, 2014. URL http://www.eumetsat.int/website/home/Satellites/ FutureSatellites/EUMETSATPolarSystemSecondGeneration/ index.html.