Satellite Data Compression
Bormin Huang Editor
Satellite Data Compression
Editor Bormin Huang Space Science and Engineering Center University of Wisconsin – Madison Madison, WI, USA [emailprotected]
ISBN 978-1-4614-1182-6 e-ISBN 978-1-4614-1183-3 DOI 10.1007/978-1-4614-1183-3 Springer New York Dordrecht Heidelberg London Library of Congress Control Number: 2011939205 # Springer Science+Business Media, LLC 2011
All rights reserved. This work may not be translated or copied in whole or in part without the written permission of the publisher (Springer Science+Business Media, LLC, 233 Spring Street, New York, NY 10013, USA), except for brief excerpts in connection with reviews or scholarly analysis. Use in connection with any form of information storage and retrieval, electronic adaptation, computer software, or by similar or dissimilar methodology now known or hereafter developed is forbidden. The use in this publication of trade names, trademarks, service marks, and similar terms, even if they are not identified as such, is not to be taken as an expression of opinion as to whether or not they are subject to proprietary rights. Printed on acid-free paper Springer is part of Springer Science+Business Media (www.springer.com)
Contents
1
2
3
4
Development of On-Board Data Compression Technology at Canadian Space Agency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Shen-En Qian
1
CNES Studies for On-Board Compression of High-Resolution Satellite Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Carole Thiebaut and Roberto Camarero
29
Low-Complexity Approaches for Lossless and Near-Lossless Hyperspectral Image Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Andrea Abrardo, Mauro Barni, Andrea Bertoli, Raoul Grimoldi, Enrico Magli, and Raffaele Vitulli FPGA Design of Listless SPIHT for Onboard Image Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Yunsong Li, Juan Song, Chengke Wu, Kai Liu, Jie Lei, and Keyan Wang
47
67
5
Outlier-Resilient Entropy Coding. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Jordi Portell, Alberto G. Villafranca, and Enrique Garcı´a-Berro
6
Quality Issues for Compression of Hyperspectral Imagery Through Spectrally Adaptive DPCM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bruno Aiazzi, Luciano Alparone, and Stefano Baronti
115
Ultraspectral Sounder Data Compression by the Prediction-Based Lower Triangular Transform . . . . . . . . . . . . . . . . . . . . . Shih-Chieh Wei and Bormin Huang
149
7
8
Lookup-Table Based Hyperspectral Data Compression . . . . . . . . . . . . Jarno Mielikainen
87
169
v
vi
9
10
11
12
13
14
Contents
Multiplierless Reversible Integer TDLT/KLT for Lossy-to-Lossless Hyperspectral Image Compression. . . . . . . . . . . Jiaji Wu, Lei Wang, Yong Fang, and L.C. Jiao
185
Divide-and-Conquer Decorrelation for Hyperspectral Data Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ian Blanes, Joan Serra-Sagrista`, and Peter Schelkens
215
Hyperspectral Image Compression Using Segmented Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Wei Zhu, Qian Du, and James E. Fowler
233
Fast Precomputed Vector Quantization with Optimal Bit Allocation for Lossless Compression of Ultraspectral Sounder Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bormin Huang
253
Effects of Lossy Compression on Hyperspectral Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chulhee Lee, Sangwook Lee, and Jonghwa Lee
269
Projection Pursuit-Based Dimensionality Reduction for Hyperspectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Haleh Safavi, Chein-I Chang, and Antonio J. Plaza
287
Contributors
Andrea Abrardo Dip. di Ingegneria dell’Informazione, Universita` di Siena, Siena, Italy Bruno Aiazzi IFAC-CNR, Via Madonna del Piano 10, 50019 Sesto F.no, Italy Luciano Alparone Department of Electronics & Telecommunications, University of Florence, Via Santa Marta 3, 50139 Florence, Italy Mauro Barni Dip. di Ingegneria dell’Informazione, Universita` di Siena, Siena, Italy Stefano Baronti IFAC-CNR, Via Madonna del Piano 10, 50019 Sesto F.no, Italy Andrea Bertoli Carlo Gavazzi Space S.p.A., Milan, Italy Ian Blanes Universitat Auto`noma de Barcelona, E-08290 Cerdanyola del Valle`s (Barcelona), Spain Roberto Camarero Ctr. National d’E´tudes Spatiales (CNES), Toulouse, France Chein-I Chang Remote Sensing Signal and Image Processing Laboratory, Department of Computer Science and Electrical Engineering, University of Maryland, Baltimore, MD, USA Department of Electrical Engineering, National Chung Hsing University, Taichung, Taiwan Qian Du Department of Electrical and Computer Engineering, Mississippi State University, Starkville, USA Yong Fang College of Information Engineering, Northwest A&F University, yangling, China James E. Fowler Department of Electrical and Computer Engineering, Mississippi State University, Starkville, USA
vii
viii
Contributors
Enrique Garcı´a-Berro Institut d’Estudis Espacials de Catalunya, Barcelona, Spain Departament de Fı´sica Aplicada, Universitat Polite`cnica de Catalunya, Castelldefels, Spain Raoul Grimoldi Carlo Gavazzi Space S.p.A., Milan, Italy Bormin Huang Space Science and Engineering Center, University of Wisconsin–Madison, Madison, WI, USA L.C. Jiao School of Electronic Engineering, Xidian University, xi’an, China Chulhee Lee Yonsei University, South Korea Jonghwa Lee Yonsei University, South Korea Sangwook Lee Yonsei University, South Korea Jie Lei State Key Laboratory of Integrated Service Networks, Xidian University, Xi’an, China Yunsong Li State Key Laboratory of Integrated Service Networks, Xidian University, Xi’an, China Kai Liu State Key Laboratory of Integrated Service Networks, Xidian University, Xi’an, China Enrico Magli Dip. di Elettronica, Politecnico di Torino, Torino, Italy Jarno Mielikainen School of Electrical and Electronic Engineering, Yonsei University, Seoul, South Korea Antonio J. Plaza Department of Technology of Computers and Communications, University of Extremadura, Escuela Politecnica de Caceres, Caceres, SPAIN Jordi Portell Departament d’Astronomia i Meteorologia/ICCUB, Universitat de Barcelona, Barcelona, Spain Institut d’Estudis Espacials de Catalunya, Barcelona, Spain Shen-En Qian Canadian Space Agency, St-Hubert, QC, Canada Haleh Safavi Remote Sensing Signal and Image Processing Laboratory, Department of Computer Science and Electrical Engineering, University of Maryland, Baltimore, MD, USA Peter Schelkens Department of Electronics and Informatics, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium Interdisciplinary Institute for Broadband Technology, Gaston Crommenlaan 8, b102, B-9050 Ghent, Belgium
Contributors
Joan Serra-Sagrista` Universitat Auto`noma de Barcelona, E-08290 Cerdanyola del Valle`s (Barcelona), Spain Juan Song State Key Laboratory of Integrated Service Networks, Xidian University, Xi’an, China Carole Thiebaut Ctr. National d’E´tudes Spatiales (CNES), Toulouse, France Alberto G. Villafranca Institut d’Estudis Espacials de Catalunya, Barcelona, Spain Departament de Fı´sica Aplicada, Universitat Polite`cnica de Catalunya, Castelldefels, Spain Raffaele Vitulli European Space Agency – ESTEC TEC/EDP, Noordwijk, The Netherlands Keyan Wang State Key Laboratory of Integrated Service Networks, Xidian University, Xi’an, China Lei Wang School of Electronic Engineering, Xidian University, xi’an, China Shih-Chieh Wei Department of Information Management, Tamkang University, Tamsui, Taiwan Chengke Wu State Key Laboratory of Integrated Service Networks, Xidian University, Xi’an, China Jiaji Wu School of Electronic Engineering, Xidian University, xi’an, China Wei Zhu Department of Electrical and Computer Engineering, Mississippi State University, Starkville, USA
ix
Chapter 1
Development of On-Board Data Compression Technology at Canadian Space Agency Shen-En Qian
Abstract This chapter reviews and summarizes the researches and developments on data compression techniques for satellite sensor data at the Canadian Space Agency in collaboration with its partners in other government departments, academia and Canadian industry. This chapter describes the subject matters in the order of the following sections.
1 Review of R&D of Satellite Data Compression at the Canadian Space Agency The Canadian Space Agency (CSA) began developing data compression algorithms as an enabling technology for a hyperspectral satellite in the 1990s. Both lossless and lossy data compression techniques have been studied. Compression techniques for operational use have been developed [1–16]. The focus of the development is the vector quantization (VQ) based near lossless data compression techniques. Vector quantization is an efficient coding technique for data compression of hyperspectral imagery because of its simplicity and its vector nature for preservation of the spectral signatures associated with individual ground samples in the scene. Reasonably high compression ratios (>10:1) can be achieved with significantly high compression fidelity for most of remote sensing applications and algorithms.
S.-E. Qian (*) Canadian Space Agency, St-Hubert, QC, Canada e-mail: [emailprotected] B. Huang (ed.), Satellite Data Compression, DOI 10.1007/978-1-4614-1183-3_1, # Springer Science+Business Media, LLC 2011
1
2
1.1
S.-E. Qian
Lossless Compression
In the 1990s, a prediction-based approach was adopted. It used a linear or nonlinear predictor in spatial or/and spectral domain to generate prediction and applies DPCM to generate residuals followed by an entropy coding. A total of 99 fixed coefficient predictors covering 1D, 2D and 3D versions were investigated [15]. An adaptive predictor that chose the best predictor from a pre-selected predictor bank was also investigated. The Consultative Committee for Space Data System (CCSDS) recommended lossless algorithm (mainly the Rice algorithm) [17] was selected as the entropy encoder to encode the prediction residuals. In order to evaluate the performance of the CCSDS lossless algorithm, an entropy encoder referred to as Base-bit Plus Overflow-bit Coding (BPOC) [18] was also utilized to compare the entropy coding efficiency. Three hyperspectral datacubes acquired using both the Airborne Visible/Near Infrared Imaging Spectrometer (AVIRIS) and the Compact Airborne Spectrographic Imager (casi) were tested. Two 3D predictors, which use five and seven nearest neighbor pixels in 3D, produced the best residuals. Compression ratios around 2.4 were achieved after encoding the residuals produced by using these two entropy encoders. All other predictors produced compression ratios smaller than 2:1. The BPOC slightly outperformed the CCSDS lossless algorithm.
1.2
Wavelet Transform Lossy Data Compression
We have also developed a compression technique for hyperspectral data using wavelet transform techniques [16]. The method for creating zerotrees as proposed by Shapiro [19] has been modified for use with hyperspectral data. An optimized multi-level lookup table was introduced which improves the performance of an embedded zerotree wavelet algorithm. In order to evaluate the performance of the algorithm, this algorithm was compared to SPIHT [20] and JPEG. The new algorithm was found to perform as well as or to surpass published algorithms which are much more computationally complex and not suitable for compression of hyperspectral data. As in Sect. 1.1 testing was done with both AVIRIS and casi data and compression ratios over 32:1 were obtained with fidelities greater than 40.0 dB.
1.3
Vector Quantization Data Compression
We selected VQ technique for data compression of hyperspectral imagery because of its simplicity and because its vector nature can preserve the spectra associated with individual ground samples in the scene. Vector quantization is an efficient coding technique. The generalized Lloyed algorithm (GLA) (sometimes also called LBG) is the most widely used VQ algorithm for image compression [21].
1 Development of On-Board Data Compression Technology. . .
3
In our application, we associate the vectors used in the VQ algorithm with the full spectral dimension of each ground sample of a hyperspectral datacube that we refer to as a spectral vector. Since the number of targets in the scene of a datacube is limited, the number of trained spectral vectors (i.e. codevectors) can be much smaller than the total number of spectral vectors of the datacube. Thus, we can represent all the spectral vectors using a codebook with comparatively few codevectors and achieve good reconstruction fidelity. VQ compression techniques make good use of the high correlation often found between bands in the spectral domain and achieves a high compression ratio. However, a big challenge to VQ compression techniques for hyperspectral imagery in terms of operational use is that it requires large computational resources, particularly for the codebook generation phase. Since the size of hyperspectral datacubes can be hundreds of times larger than those for traditional remote sensing, the processing time required to train a codebook or to encode a datacube using the codebook could also be tens to hundreds of times larger. In some applications, the problem of training time can largely be avoided by training a codebook only once, and henceforth applying it repeatedly to all subsequent datacubes to be compressed as adopted in the conventional 2D image compression. This works well when the datacube to be compressed is bounded by the training set used to train the codebook. However, in hyperspectral remote sensing, it is in general very difficult to obtain a so-called “universal” codebook that spans many datacubes to the required degree of fidelity. This is partly because the characteristics of the targets (season, location, illumination, view angle, the needs for accurate atmospheric effects) and instrument configuration (spectral and spatial resolution, spectral range, SNR of the instruments) introduce high variability in the datacubes and partly because of the need for high reconstruction fidelity in its downstream use. For these reasons, it is preferred that a new codebook is generated for every datacube that is to be compressed, and is transmitted to the decoder together with the index map as the compressed data. Thus, the main goal in the development of VQ based techniques to compress hyperspectral imagery is to seek a much faster and more efficient compression algorithm to overcome this challenge, particularly for on-board applications. In this chapter, the compression of a hyperspectral datacube using the conventional GLA algorithm with a new trained codebook for a datacube is referred to as 3DVQ herein. The CSA developed an efficient method of representing spectral vectors of hyperspectral datacubes referred to as Spectral Feature Based Binary Code (SFBBC) [1]. With SFBBC code, Hamming distance, rather than Euclidean distance, could be used in codebook training and codevector matching (coding) in the process of the VQ compression for hyperspectral imagery. The Hamming distance is a simple sum of logical bit-wise exclusive-or operations. It is much faster than Euclidean distance. The SFBBC based VQ compression technique can speed up the compression process by a factor of 30–40 at a fidelity loss of PSNR . Equivalently, each pixel value g(n) can be predicted as a fuzzy switching, i.e., a blending, of the outputs of all the predictors, which are defined as ^ gm ðnÞ ¼ with the fuzzy prediction, g^(n), given by " # M X ~ m ðnÞ ^gm ðnÞ g^ðnÞ ¼ round U
(6.5)
(6.6)
m¼1
The right term of (6.6) is rounded to integer to yield integer valued prediction errors, i.e., eðnÞ ¼ gðnÞ g^ðnÞ, that are sent to the entropy coding section. Concerning S-RLP, once blocks have been classified and labeled, together with the attached optimized predictor, each band is raster scanned and predictors are activated based on the classes of crossed blocks. Thus, each pixel belonging to one block of the original partition, g(x, y), is predicted as a g^(x, y) by using the one out of the M predictors that was found to better fit the statistics of that class of data block in the MMSE sense. The integer valued prediction errors, viz., eðx; yÞ ¼ gðx; yÞ round½^ gðx; yÞ, are delivered to the context-coding section, identical to that of FMP.
4.4
Context Based Arithmetic Coding
As evidenced in Fig. 6.2, a more efficient coding of the prediction errors can be obtained if a context classification of the residuals is put ahead the arithmetic coding block. A context function may be defined and measured on prediction errors laying within a causal neighborhood, possibly larger than the prediction support, as the RMS value of prediction errors (RMSPE). The context function should capture the non-stationarity of prediction errors, regardless of their spatial correlation. Again, causality of neighborhood is necessary in order to make the same information available both at the encoder and at the decoder. At the former, the probability density function (PDF) of RMSPE is calculated and partitioned into a number of intervals chosen as equally populated; thus, contexts are equiprobable as well. This choice is motivated by the use of adaptive arithmetic coding for encoding the errors belonging to each class. Adaptive entropy coding, in general, does not require previous knowledge of the statistics of the source, but benefits from a number of data large enough for training, which happens simultaneously with coding. The source given by each class is further split into sign bit and magnitude. The former is strictly random and is coded as it stands, the latter exhibits a reduced variance in each class; thus, it may be coded with fewer bits than the original residue. It is noteworthy that such a context-coding procedure is independent of the particular
6 Quality Issues for Compression of Hyperspectral Imagery. . .
127
method used to decorrelate the data. Unlike other schemes, e.g., CALIC [52], in which context-modeling is embedded in the decorrelation procedure, the method [7] can be applied to any DPCM scheme, either lossless or near-lossless, without adjustments in the near-lossless case, as it happens to other methods [51].
5 Lossless Hyperspectral Image Compression Based on LUTs In a recent published paper [32], Mielikainen introduced a very simple prediction method, referred as LUT-NN, which predicts the current pixel by taking the value of its nearest neighbor (NN), defined as the pixel previously transmitted in the current band, having the following two properties: (1) the pixel in the previous band, at the same position of the NN pixel, has the same value of the pixel in the previous band at the same position of the current pixel, and (2) among all the pixels fulfilling the previous property, the NN pixel is the spatially closest to the current pixel along the scan path. Such an algorithm can be effectively implemented by means of a dynamically updated lookup table (LUT). The prediction value taken from the LUT is that of the cell indexed by the value of the pixel at the current pixel position in the previous band and corresponds to the value of the NN pixel, which has been inserted in this cell in a previous updating pass [32]. Surprisingly, such a simple method, notwithstanding only one previous band is considered, performs comparably to the most advanced DPCM methods on the 1997 AVIRIS data set (16-bit radiance format). However, this improvement occurs only if some particular calibrated sequences are considered. Actually, in the case of the 1997 AVIRIS data set, the calibration procedures produces for each band radiance histograms that are very unbalanced, i.e., some values are much more frequent than others. So, the LUT-NN method forces its predictions to yield the most frequent values in the band to be compressed, unlike conventional prediction strategies usually do. In any case, this artificial efficiency has suggested the investigation of more sophisticated versions of LUT-NN. A first attempt was the method proposed by Huang et al., which was named LAIS-LUT [23]. This algorithm utilizes two LUTs, respectively containing the values of the two closest NNs to the current pixel. The prediction of the current pixel is obtained by choosing the value that is more similar to a reference prediction, which takes into account the interband gain between the current and the previous bands, as it is indicated by its acronym Locally Average Interband Scaling (LAIS). The LAIS-LUT algorithm yields significant coding benefits over LUT-NN at a moderate extra cost, thanks to its better accounting of the spatial correlation, even if the exploitation of the spectral correlation is still limited to the previous band only. The idea underlying LUT-NN and LAIS-LUT can be further generalized by considering a more advanced reference prediction and by exploiting also the spectral correlation of the sequences. The adoption of the S-RLP and S-FMP methods as reference predictors has brought to the two generalizations proposed in [10] and denoted as S-RLP-LUT and S-FMP-LUT.
128
B. Aiazzi et al.
These two algorithms feature a complete exploitation of spatial and spectral correlations, because the prediction value is obtained by considering more than two LUTs for bands, say M, where M is usually chosen equal to 4, and an arbitrary number of previous bands, say N, where N may reach 20. The decision among the M N possible prediction values is based on the closeness of the candidate predictions to the S-FMP or S-RLP reference predictions, which can span the same N previous bands. The advantage of considering a large number of LUTs is strictly connected to the utilization of S-FMP and S-RLP as reference predictors. In fact, by adopting the simpler LAIS predictor as reference, the advantages of more than one band for LUT-based prediction and of more than two LUTs for each band quickly vanish as M and N increase. Figure 6.4 shows how S-FMP-LUT and S-RLP-LUT work, i.e., how an algorithm based on multiple LUTs can be connected with the S-RLP and S-FMP advanced DPCM predictors. By using the same notation of Sect. 4.3.1, let g(x, y, z) (0 g (x, y, z) < gfs) denote the intensity level of the sequence at row x (1 x Nx), column y (1 y Ny), and band z (1 z Nz). Let also Lm;n;z ½; m ¼ 1; M; n ¼ 1; N indicate the set of N M LUTs used for prediction of band z. All LUTs are of length gfs. Eventually, let g^r ðx; y; zÞ be the reference prediction for pixel (x, y) in band z, obtained by means of S-RLP or S-FMP. The multiband LUT prediction is given by g^ðx; y; zÞ ¼ g^m;^ ^ n ðx; y; zÞ
(6.7)
in which ^ n^g¼ arg min fm;
m¼1;M n¼1;N
n g^
m;n ðx; y; zÞ
o g^r ðx; y; zÞ
(6.8)
and g^m;n ðx; y; zÞ ¼ Lm;n;z ½gðx; y; z
nÞ;
m ¼ 1; M; n ¼ 1; N
(6.9)
are the N M possible prediction values among which the final prediction g^ðx; y; zÞ is chosen. The values of g^m;n ðx; y; zÞ belong to the current band and have been stored in the set of multiband LUTs during the previous updating steps. After the final prediction value has been produced according to (6.7), the set of multiple LUTs is updated, analogously to [23, 32] in the following way: Lm;n;z ½gðx; y; z nÞ ¼ Lm 1;n;z ½gðx; y; z m ¼ M; M 1; 2; n ¼ 1; N L1;n;z ½gðx; y; z nÞ ¼ gðx; y; zÞ; n ¼ 1; N:
nÞ;
(6.10)
6 Quality Issues for Compression of Hyperspectral Imagery. . .
129
a CRISP OR FUZZY ADAPTIVE PREDICTION (reference prediction) Previous Bands
FINAL LUT PREDICTION (the most fitting with reference prediction)
Current Band LUT PREDICTIONS (4 per previous band max 20 bands)
− + +
CODING OF RESIDUALS
b CRISP OR FUZZY ADAPTIVE PREDICTION (reference prediction) Previous Bands
FINAL LUT PREDICTION (the most fitting with reference prediction)
Reconstructed Current Band LUT PREDICTIONS (4 per previous band max 20 bands)
+ Z-1
+
+
DECODING OF RESIDUALS
Fig. 6.4 General scheme for the adaptation of the basic LUT method to S-RLP and S-FMP algorithms (a) encoder; (b) decoder. The crisp or fuzzy adaptive DPCM algorithms work as advanced reference predictors for an optimal selection among the multiple LUT-based prediction values
All LUTs are initialized with the gfs value that is outside the range of the data. Whenever such a value is returned by (6.7), i.e., no pixel exists fulfilling the two NN properties, the reference prediction g^r ðx; y; zÞ is used instead of g^m;^ ^ n ðx; y; zÞ. At the decoder, the set of multiband LUTs and the reference prediction are calculated from the previously decoded lossless data, by following the same procedure as at the encoder. The performances of S-FMP-LUT and S-RLP-LUT on the 1997 AVIRIS data set are impressive and shows an average improvement with respect to LAIS-LUT of
130
B. Aiazzi et al.
more than 20%. However, in the same paper [10], the authors note that LUT-based methods lose their effectiveness on unprocessed, i.e., non-calibrated, images, which do not contain calibration-induced artifacts. This limitation of the LUT-based algorithms has been strengthened in [27], where it is shown that LUT-based methods are ineffective on the 2006 AVIRIS data set. In fact, a different calibration procedure is applied to the new data, that exhibit an improved radiometric resolution. Such procedure does not generate recurrent pixel values in the radiance histogram. Finally, LUT-based methods are ineffective also in lossy compression, and in particular in near-lossless and virtually-lossless cases, because the quantization of the predictor errors smoothes the data histograms from which these algorithms take the prediction values.
6 Near-Lossless Compression So far, quantization in the FMP and RLP schemes was not addressed; that is, lossless compression was described. Quantization is necessarily introduced to allow a reduction in the code rate to be achieved [24] at the price of some distortions in the decompressed data. Although trellis coded quantization may be optimally coupled with a DPCM scheme [25], its complexity grows with the number of output levels and especially with the complexity of predictors. Therefore, in the following only linear and logarithmic quantization will be concerned. The latter is used to yield relative error bounded compression.
6.1
Distortion Measurements
Before discussing quantization in DPCM schemes, let us review the most widely used distortion measurements suitable for single-band image data (2D) and multiband image data (3D), either multispectral or hyperspectral.
6.1.1
Radiometric Distortion
Let 0 g(x, y) gfs denote an N-pixel digital image and let g~ðx; yÞ be its possibly distorted version achieved by compressing g(x, y) and decompressing the outcome bit stream. Widely used distortion measurements are reported in the following. Mean absolute error (MAE), or L1 norm, MAE ¼
1 XX jgðx; yÞ N x y
g~ðx; yÞj;
(6.11)
6 Quality Issues for Compression of Hyperspectral Imagery. . .
131
Mean Squared Error (MSE), or L22, MSE ¼
1 XX ½gðx; yÞ N x y
g~ðx; yÞ2 ;
(6.12)
Root MSE (RMSE), or L2, RMSE ¼ Signal to Noise Ratio (SNR)
pffiffiffiffiffiffiffiffiffiffi MSE;
SNRðdBÞ ¼ 10 log10
(6.13)
g2 ; 1 MSE þ 12
(6.14)
Peak SNR (PSNR) PSNRðdBÞ ¼ 10 log10
g2fs 1 MSE þ 12
;
(6.15)
Maximum absolute distortion (MAD), or peak error, or L1, MAD¼ maxx;y fjgðx; yÞ
~ gð x; yÞjg;
(6.16)
Percentage maximum absolute distortion (PMAD) PMAD¼ maxx;y
jgðx; yÞ
~ gð x; yÞj
gðx; yÞ
100:
(6.17)
In both (6.14) and (6.15) the MSE is incremented by the variance of the integer roundoff error, to handle the limit lossless case, when MSE¼0. Thus, SNR and PSNR will be upper bounded by 10 log10 ð12 g2 Þ and 10 log10(12 gfs2), respectively. When multiband data are concerned, let vz ≜ gz(x, y), z ¼ 1, , Nz, denote the zth component of the original multispectral pixel vector v and ~vz ¼D ~gz ðx; yÞ; z ¼ 1; ; N z its distorted version. Some of the radiometric distortion measurements (6.11)–(6.17) may be extended to vector data. Average vector RMSE (VRMSE), or L1(L2) (the innermost norm L2 refers to vector space (z), the outer one to pixel space (x, y)) VRMSEavg ¼
ffi X 1 X rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½gz ðx; yÞ ~gz ðx; yÞ2 ; N x;y z
(6.18)
132
B. Aiazzi et al.
Peak VRMSE, or L1(L2), VRMSEmax ¼ maxx;y
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rX ½gz ðx; yÞ ~gz ðx; yÞ2 ;
(6.19)
z
SNR SNR ¼ 10 log10 P
PSNR
x;y;z
PSNR ¼ 10 log10 P
P
x;y;z
g2z ðx; yÞ
½gz ðx; yÞ
~gz ðx; yÞ2
N N z g2fs
x;y;z
½gz ðx; yÞ
(6.20)
;
~gz ðx; yÞ2
;
(6.21)
MAD, or L1(L1),
MAD¼ maxfjgz ðx; yÞ x;y;z
~ gz ðx; yÞjg;
(6.22)
PMAD PMAD¼ max x;y;z
jgz ðx; yÞ ~ gz ðx; yÞj 100: gz ðx; yÞ
(6.23)
In practice, (6.18) and (6.19) are respectively the average and maximum of the Euclidean norm of the distortion vector. SNR (6.20) is the extension of (6.14) to the 3D data cube. PSNR is the maximum SNR, given the full-scales of each vector component. MAD (6.22) is the maximum over the pixel set of the maximum absolute component of the distortion vector. PMAD (6.23) is the maximum percentage error over each vector component of the data cube.
6.1.2
Spectral Distortion
Given two spectral vectors v and ~ v both having L components, let v ¼ {v1, v2, , vL} be the original spectral pixel vector vz ¼ gz(x, y) and ~v ¼ f~v1 ; ~v2 ; ; ~vL g its distorted version obtained after lossy compression and decompression, i.e., ~vz ¼ ~ gz ðx; yÞ. Analogously to the radiometric distortion measurements, spectral distortion measurements may be defined. The spectral angle mapper (SAM) denotes the absolute value of the spectral angle between the couple of vectors: D arccos SAMðv; ~ vÞ ¼ (6.24) jjvjj2 jj~vjj2
6 Quality Issues for Compression of Hyperspectral Imagery. . .
133
in which < , > stands for scalar product. SAM can be measured in either degrees or radians. Another measurement especially suitable for hyperspectral data (i.e., for data with large number of components) is the spectral information divergence (SID) [20] derived from information-theoretic concepts: SIDðv; ~ vÞ ¼ Dðvjj~ vÞ þ Dð~vjjvÞ
(6.25)
with Dðvjj~ vÞ being the Kullback-Leibler distance (KLD), or entropy divergence, or discrimination [24], defined as D Dðvjj~ vÞ ¼
L X z¼1
p pz log z qz
(6.26)
in which D pz ¼
vz jjvjj1
vz D ~ and qz ¼ jj~vjj
(6.27) 1
In practice SID is equal to the symmetric KLD and can be compactly written as SIDðv; ~ vÞ ¼
L X z¼1
ðpz
p qz Þ log z qz
(6.28)
which turns out to be symmetric, as one can easily verify. It can be proven as well that SID is always nonnegative, being zero iff. pz qz, 8z, i.e., if v is parallel to ~v. The measure unit of SID depends on the base of logarithm: nat/vector with natural logarithms and bit/vector with logarithms in base two. Both SAM (6.24) and SID (6.28) may be either averaged on pixel vectors, or the maximum may be taken instead, as more representative of spectral quality. Note that radiometric distortion does not necessarily imply spectral distortion. Conversely, spectral distortion is always accompanied by a radiometric distortion, that is minimal when the couple of vectors have either the same Euclidean length (L2) for SAM, or the same city-block length (L1), for SID.
6.2
Linear and Logarithmic Quantization
In this subsection, linear and logarithmic quantization are described. Quantization has the objective of reducing the transmission rate. Linear quantization is capable of controlling near-lossless coding while logarithmic quantization is used to yield relative error bounded compression.
134
6.2.1
B. Aiazzi et al.
Linear Quantization
In order to achieve reduction in bit rate within the constraint of a near-lossless compression, prediction errors are quantized, with a quantization noise feedback loop embedded into the encoder, so that the current pixel prediction is formulated from the same “noisy” data that will be available at the decoder, as shown in Fig. 6.1a,b. Prediction errors, e(n) ≜ g(n) g^ðnÞ, may be linearly quantized with a step size D as eD ðnÞ ¼ round½eðnÞ=D and delivered to the context-coding section, as shown in Fig. 6.1a. The operation of inverse quantization, e~ðx; yÞ ¼ eD ðx; yÞ D introduces an error, whose variance and maximum modulus are D2 / 12 and bD / 2c, respectively. Since the MSE distortion is a quadratic function of the D, an odd-valued step size yields a lower L1 distortion for a given L2 than an even size does; thus, odd step sizes are preferred for near-lossless compression. The relationship between the target peak error, i.e., e 2 Zþ , and the step size to be used is D ¼ 2e þ 1.
6.2.2
Logarithmic Quantization
For the case of a relative-error bounded compression a rational version of prediction error must be envisaged. Let us define the relative prediction error (RPE) as ratio of original to predicted pixel value: D gðnÞ rðnÞ¼ g^ðnÞ
(6.29)
The rational nature of RPE, however, makes linear quantization unable to guarantee a strictly user-defined relative-error bounded performance. Given a step size D 2 R (D > 0, D6¼1), define direct and inverse logarithmic quantization (Log-Q) of t 2 R, t > 0, as D round½log ðtÞ ¼ round½logðtÞ= logðDÞ QD ðtÞ¼ D
QD 1 ðlÞ ¼ Dl
(6.30)
Applying (6.30) to (6.29) yields QD ½rðnÞ ¼ round
logðgðnÞÞ logð^ gðnÞÞ log D
(6.31)
Hence, Log-Q of RPE is identical to Lin-Q of logðgðnÞÞ logð^ gðnÞÞ with a step size logD. If a Log-Q with a step size D is used to encode pixel RPE’s (6.29), it can be proven that the ratio of original to decoded pixel value is strictly bounded around one
6 Quality Issues for Compression of Hyperspectral Imagery. . .
pffiffiffiffi 1 min D; pffiffiffiffi D
pffiffiffiffi 1 gðnÞ max D; pffiffiffiffi g~ðnÞ D
135
(6.32)
From (6.32) and (6.29) it stems that the percentage pixel distortion is upper bounded PMAD ¼ max
pffiffiffiffi D
1; 1
1 pffiffiffiffi 100 D
(6.33)
depending on whether D > 1, or 0 < D < 1. Hence, the relationship between the target percentage peak error, r, and the step size will be, e.g., for D > 0, D ¼ ð1 þ r=100Þ2 .
7 Virtually Lossless Compression The term virtually lossless indicates that the distortion introduced by compression appears as an additional amount of noise, besides the intrinsic observation noise, being statistically independent of the latter, as well as of the underlying signal. Its first order distribution should be such that the overall probability density function (PDF) of the noise corrupting the decompressed data, i.e., intrinsic noise plus compression-induced noise, closely matches the noise PDF of the original data. This requirement is trivially fulfilled if compression is lossless, but may also hold if the difference between uncompressed and decompressed data exhibits a peaked and narrow PDF without tails, as it happens for near lossless techniques, whenever the user defined MAD is sufficiently smaller than the standard deviation sn of the background noise. Both MAD and sn are intended to be expressed in either physical units, for calibrated data, or as digital counts otherwise. Therefore, noise modeling and estimation from the uncompressed data becomes a major task to accomplish a virtually lossless compression [6]. The underlying assumption is that the dependence of the noise on the signal is null, or weak. However, signal independence of the noise may not strictly hold for hyperspectral images, especially for newgeneration instruments. This further uncertainty in the noise model may be overcome by imposing a margin on the relationship between target MAD and RMS value of background noise. For a DPCM scheme, the relationship between MAD and quantization step size is D ¼ 2MAD + 1, while the relationship between the variance of quantization noise, which is equal to MSE, and the step size is e 2 ¼ D2 =12. Hence, the rationale of virtually lossless compression can be summarized as follows. Firstly, measure the background noise RMS, sn; if sn < 1, lossless compression is mandatory. If 1 sn < 3, lossless compression is recommended, but near lossless compression with MAD ¼ 1 (D ¼ 3) is feasible. For 3 sn < 5, a strictly virtually lossless compression would require MAD ¼ 1, and so on.
136
B. Aiazzi et al.
The signal may have been previously quantized based on different requirement; afterwards a check on the noise is made to decide whether lossless compression is really necessary, or near lossless compression could be used instead without penalty, being “de facto” virtually lossless. The key to achieve a compression preserving the scientific quality of the data for remote-sensing is represented by the following twofold recommendation: 1. Absence of tails in the PDF of the error between uncompressed and decompressed image, in order to maximize the ratio RMSE / MAD, i.e., e/MAD, or equivalently to minimize MAD for a given RMSE. 2. MSE lower by one order of magnitude (10 dB) than the variance of background noise sn2. Near-lossless methods are capable to fulfill such requirements, provided that the quantization step size D is chosen as an odd integer such that D sn. More exactly, the relationship between MAD and sn, also including a margin of approximately 1 dB, is: MAD ¼ bmaxf0; ðsn
1Þ=2gc
(6.34)
Depending on application context and type of data, the relationship (6.34) may be relaxed, e.g., by imposing that the ratio sn2 / e2 is greater than, say, 3 dB, instead of the 11 dB, given by (6.34). If the data are intrinsically little noisy, the protocol may lead to the direct use of lossless compression, i.e., D ¼ 1. If compression ratios greater than those of the reversible case are required, near lossless compression with MAD 1 of low-noise bands may be inadequate to preserve the scientific quality of the data, because the compression-induced MSE is not one order of magnitude lower than sn2, as it would be recommended for virtually lossless compression. However, it may become mandatory to increase compression ratios above the values typical of the strictly virtually lossless protocol. To adjust the compression ratio, a real valued positive scale factor q can be introduced, such that the quantization step size of the nth band is given by: Dn ¼ round½q sn :
(6.35)
where the roundoff is to the nearest odd integer. If q 1 a strictly virtually lossless compression is achieved, since the compression-induced quadratic distortion is less than one tenth of the intrinsic noisiness of the data. Otherwise, if q > 1, compression is widely virtually lossless, even though distortion is properly allocated among the spectral bands. As an example, Fig. 6.5 shows quantization step sizes for three different values of q, if the test sequence Cuprite Mine of AVIRIS 1997 is near lossless compressed by means of the S-RLP algorithm. The estimation of the background noise RMS of the test sequence has been performed by considering a method based on the joint 2D PDF of the local statistics, described in [4]. To compare virtually lossless compression with a unique quantizer for the whole data cube, the step size of the latter, yielding the same
6 Quality Issues for Compression of Hyperspectral Imagery. . .
137
q=0.5 q=1 q=2
250
200
∆
150
100
50
0 0
50
100
150
200
Band number
Fig. 6.5 Quantization step size varying with band number for the test sequence Cuprite Mine
compression ratio as the former, is the odd integer roundoff of the geometric mean of the step sizes in (6.35), at least for small distortions, in order to have an independent quantizer at each band.
8 Experimental Results and Comparisons Experimental results are presented for the two data sets constituted by AVIRIS 1997 and AVIRIS 2006 data. Lossless compression is first considered on both data sets. Afterwards some results are presented for near-lossless and virtually lossless compression on a particular scene of the Cuprite Mine sequence of the AVIRIS 1997 data set. Eventually, some considerations on virtually lossless compression on discrimination of materials are introduced.
8.1
Lossless Compression Performance Comparisons
The lossless compression experiments have been performed by comparing some of the most performing schemes reported in the literature such as FL [27], SLSQ [42], M-CALIC [30], C-DPCM [33], S-RLP and S-FMP for the classical methods, and LUT-NN [32], LAIS-LUT [23], S-RLP-LUT and S-FMP-LUT for the LUT-based algorithms. Eventually, also TSP [27] is reported, an algorithm specifically designed to cope with the calibrated-induced artifacts. Some of these methods have not been made available by authors on the 2006 AVIRIS data set.
138
B. Aiazzi et al.
Table 6.1 Bit rates (bit/pel/band on disk) for lossless compression of AVIRIS 1997 test hyperspectral images. S-RLP-LUT and S-FMP-LUT use N ¼ 20 previous bands and M ¼ 4 LUT’s per band. Best results are reported in bold Cuprite Jasper Lunar Moffett Average FL 4.82 4.87 4.83 4.93 4.86 SLSQ 4.94 4.95 4.95 4.98 4.96 M-CALIC 4.89 4.97 4.80 4.65 4.83 C-DPCM 4.68 4.62 4.75 4.62 4.67 S-RLP 4.69 4.65 4.69 4.67 4.67 S-FMP 4.66 4.63 4.66 4.63 4.64 LUT-NN 4.66 4.95 4.71 5.05 4.84 LAIS-LUT 4.47 4.68 4.53 4.76 4.71 S-RLP-LUT 3.92 4.05 3.95 4.09 4.00 S-FMP-LUT 3.89 4.03 3.92 4.05 3.97 TSP 3.77 4.08 3.81 4.12 3.95
8.1.1
1997 AVIRIS Data Set
Table 6.1 lists the bit rates for several lossless compression methods applied to the standard 1997 AVIRIS images. The compared methods are FL [27], SLSQ [42], M-CALIC [30], C-DPCM [33], S-RLP, S-FMP, LUT-NN [32], LAIS-LUT [23], S-RLP-LUT, S-FMP-LUT and TSP [27]. It is apparent from the results reported in Table 6.1 that the best average performance is obtained by TSP and by the LUT-based algorithms, which are very effective even if not specifically designed for this type of data, because they are able to exploit the artificial regularities in the image histogram. Concerning the S-RLP algorithm, Fig. 6.6 reports a comparison with 3D-RLP by varying the number S of the pixels in the prediction causal neighborhood P S ðnÞ. In this experiment, the fourth scene of the test image Cuprite ’97 has been reversibly compressed by means of S-RLP and 3D-RLP. Bit rates, reported in bits per pixel per band, include all coding overhead: predictor coefficients, block labels, and arithmetic codewords. For each l, prediction lengths equal to 5, 14, and 20 have been considered together with 16 predictors. 3D prediction of RLP is always carried out from a couple of previous bands, except for the first band, coded in intra mode, i.e., by 2D DPCM, and the second band, which is predicted from one previous band only. The purely spectral 1D prediction of S-RLP is carried out from the available previous bands up the requested prediction length. As it appears from the plots in Fig. 6.6, S-RLP outperforms 3D-RLP, especially when prediction length is lower, a case of interest for customized satellite on-board implementations. The performance of both S-RLP and 3D-RLP cannot be improved significantly by increasing the number and length of predictors, because of overhead information increasing as well.
6 Quality Issues for Compression of Hyperspectral Imagery. . .
139
5 3D−RLP S−RLP
4.95
Bit Rate (bit/pel/band)
4.9 4.85 4.8 4.75 4.7 4.65 4.6 4.55 4.5
5 pels
14 pels
20 pels
Prediction Length
Fig. 6.6 S-RLP versus 3D-RLP varying with the prediction length. Bit rates on disk for the lossless compression of the fourth scene of AVIRIS Cuprite Mine 1997
Fig. 6.7 Band 50 of the 2006 AVIRIS Yellowstone 10 scene
8.1.2
2006 AVIRIS Data Set
A new data set of calibrated and uncalibrated AVIRIS images has been provided by Consultative Committee for Space Data Systems (CCSDS) and is now available to scientific users for compression experiments at the web site of address:(http:// compression.jpl.nasa.gov/hyperspectral). This data set consists of five calibrated and the corresponding raw 16-bit images acquired over Yellowstone, WY. Each image is composed by 224 bands and each scene (the scene numbers are 0, 3, 10, 11, 18) has 512 lines [27]. The 50th band of Yellowstone 10 is reported in Fig. 6.7.
140
B. Aiazzi et al.
Table 6.2 Bit rates (bit/pel/band on disk) for lossless compression of the calibrated 2006 AVIRIS data set. S-RLP-LUT and S-FMP-LUT have been utilized only for the first two scenes. The best results are reported in bold Yellowstone Yellowstone Yellowstone Yellowstone Yellowstone 18 Average 11 10 3 0 FL 3.91 3.79 3.37 3.59 3.90 3.71 S-RLP 3.58 3.43 2.95 3.27 3.46 3.34 S-FMP 3.54 3.39 2.94 3.25 3.44 3.31 LUT-NN 4.82 4.62 3.96 4.34 4.84 4.52 LAIS-LUT 4.48 4.31 3.71 4.02 4.48 4.20 SRLP-LUT 3.95 3.82 n.a. n.a. n.a. n.a. S-FMP-LUT 3.91 3.78 n.a. n.a. n.a. n.a. TSP 3.99 3.86 3.42 3.67 3.97 3.78
Table 6.3 Bit rates (bit/pel/band on disk) for lossless compression of the uncalibrated 2006 AVIRIS data set. S-RLP-LUT and S-FMP-LUT have been utilized only for the first two scenes. The best results are reported in bold Yellowstone Yellowstone Yellowstone Yellowstone Yellowstone 18 Average 11 10 3 0 FL 6.20 6.07 5.60 5.81 6.26 5.99 S-RLP 5.88 5.72 5.21 5.54 5.75 5.62 S-FMP 5.84 5.67 5.18 5.48 5.68 5.57 LUT-NN 7.14 6.91 6.26 6.69 7.20 6.84 LAIS-LUT 6.78 6.60 6.00 6.30 6.82 6.50 SRLP-LUT 6.21 6.05 n.a. n.a. n.a. n.a. S-FMP-LUT 6.17 6.01 n.a. n.a. n.a. n.a. TSP 6.27 6.13 5.64 5.88 6.32 6.05
Tables 6.2 and 6.3 report the lossless compression performances on the 2006 AVIRIS calibrated and uncalibrated data sets, respectively. The compared methods are FL [28], S-RLP, S-FMP, LUT-NN [32], LAIS-LUT [23], TSP [27], S-RLPLUT and S-FMP-LUT. For these last algorithms, only two scenes have been compressed. The scores show that the LUT-based methods and TSP are not able to obtain the same performances on both calibrated and uncalibrated 2006 AVIRIS data set with respect to the best performances obtained on the 1997 AVIRIS data set. S-RLP and S-FMP are the most effective and gain more than 10% and 6% over FL and 12% and 7% on TSP, on calibrated and uncalibrated data, respectively. FL algorithm is effective notwithstanding its simplicity.
8.2
Near Lossless Compression Performance Comparisons
Near lossless compression tests have been performed on the fourth scene of Cuprite Mine. Rate-Distortion (RD) plots are reported in Fig. 6.8a for S-RLP and 3D-RLP operating with M ¼ 16 predictors and S ¼ 20 coefficients per predictor. PSNR of
6 Quality Issues for Compression of Hyperspectral Imagery. . .
a
b 105
16
100
14
10
MAD
90 85
8 6
3D-RLP SRLP
80
4
75 70
3D-RLP SRLP
12
95
PSNR (dB)
141
2 0 0
1
2
3
4
Bit Rate (bit/pel per band)
5
1
2
3
4
5
Bit Rate (bit/pel per band)
Fig. 6.8 Performance plots for near-lossless compression of AVIRIS CupriteMine ’97 test hyperspectral images: (a) PSNR vs. bit rate; (b) maximum absolute distortion (MAD) vs. bit rate
the whole image is calculated from the average MSE of the sequence of bands. Due to the sign bit, the full scale gfs in (6.15) was set equal to 215 1 ¼ 32,767 instead of 65,535, since small negative values, introduced by removal or dark current during calibration, are very rare and totally missing in some scenes. Hence, the PSNR attains a value of 10log10(12gfs2) 102 dB, due to integer roundoff noise only, when reversibility is reached. The correction for roundoff noise has a twofold advantage. Firstly, lossless points appear inside the plot and can be directly compared. Secondly, all PSNR-bit rate plots are straight lines with slope 6 dB/bit for bit rates larger than, say, 1 bit/pel, in agreement with RD theory [24] (with a uniform threshold quantizer). For lower bit rates, the quantization noise feedback causes an exponential drift from the theoretical straight line. The results follow the same trends as the lossless case for S-RLP and are analogous to those of 3D-RLP, reported in [6], and of M-CALIC, reported in [30]. The near-lossless bit rate profiles are rigidly shifted downward from the lossless case by amounts proportional to the logarithm of the quantization-induced distortion. This behavior does not occur for low bit rates, because of the quantization noise feedback effect: prediction becomes poorer and poorer as it is obtained from the highly distorted reconstructed samples used by the predictor, which must be aligned to the decoder. Interestingly, the difference in bit rate between S-RLP and 3D-RLP at a given PSNR is only two hundredths of bit/pel near the lossless point, but grows up to one tenth of bit/pel at a rate equal to 1.5 bit/pel, typical of a high-quality lossy compression. Comparisons with an up-to-date algorithm [35], implementing the state-of-the-art JPEG2000 multi-component approach, reveal that S-RLP outperforms the wavelet-based encoder by approximately 3.5 dB at 2 bit/pel. However, this difference reduces to 2.5 dB at 1 bit/pel and vanishes around 0.25 bit/pel, because of the quantization noise feedback effect, which is missing in the 3D wavelet coder. This moderate loss of performance is the price that
142
B. Aiazzi et al.
a
b 12
12
Bit rate (bit/pixel)
10 8 6 4 2 0
PMAD = 0.075 PMAD = 0.125 PMAD = 0.1875 PMAD = 0.4
10
Bit rate (bit/pixel)
MAD = 1 MAD = 3 MAD = 5 MAD = 7 MAD = 9 MAD = 11
8 6 4 2
50
100
150
Spectral band number
200
50
100
150
200
Spectral band number
Fig. 6.9 Bit rates produced by 3D RLP on the Cuprite sequence of bands: (a) linear quantization to yield user-defined MAD values; (b) logarithmic quantization to yield user-defined PMAD values
embedded coders have to pay. DPCM does not allow progressive reconstruction to be achieved, but yields higher PSNR, at least for medium-high bit rates. The further advantage of DPCM is that it is near-lossless, unlike JPEG 2000, which can be made lossless, but not near-lossless, unless an extremely cumbersome quantizer is employed, with further loss in performances [15]. The near-lossless performance is shown in the MAD-bit rate plots of Fig. 6.8b. Since the average standard deviation of the noise was found to be around 10 according to [3], a virtually-lossless compression (maximum compression-induced absolute error lower than the standard deviation of the noise [29]) is given by S-RLP at a bit rate around 1.6 bit/pel/band, thereby yielding a compression ratio CR ¼ 10 with respect to uncompressed data and CR 3 relative to lossless compression. Figure 6.8 evidences that the increment in performance of S-RLP over 3D-RLP is more relevant for such a bit rate than for higher rates. A different experiment reported for 3D-RLP but representative of all DPCMbased methods is reported in Fig. 6.9a for linear quantization and in Fig. 6.9b for the logarithmic quantizer. All bands have been compressed in both MAD-constrained mode (linear quantization) and PMAD constrained mode (logarithmic quantization). Bit rates varying with band number, together with the related distortion parameters are shown in Fig. 6.9. The bit rate plots follow similar trends varying with the amount of distortion, but quite different trends for the two types of distortion (i.e., either MAD or PMAD). For example, around the water vapor absorption wavelengths ( Band 80) the MAD-bounded plots exhibit pronounced valleys, that can be explained because the intrinsic SNR of the data becomes lower; thus the linear quantizer dramatically abates the noisy prediction errors. On the other hand, the PMAD-bounded encoder tends to quantize the noisy residuals more finely when the signal is lower. Therefore bit rate peaks are generated instead of valleys. More generally speaking, bit rate peaks from the PMAD-bounded encoder are
6 Quality Issues for Compression of Hyperspectral Imagery. . .
a
b 300
0.18 PMAD−bounded (max) PMAD−bounded (avg.) MAD−bounded (max) MAD−bounded (avg.)
PMAD−bounded (max) PMAD−bounded (avg.) MAD−bounded (max) MAD−bounded (avg.)
0.16 0.14
SAM (degrees)
250 200
V−RMSE
143
150 100
0.12 0.1 0.08 0.06 0.04
50
0.02 0 1
1.5
2
2.5
3
3.5
4
Average bit rate (bit per pixel per band)
4.5
0 1
1.5
2
2.5
3
3.5
4
4.5
Average bit rate (bit per pixel per band)
Fig. 6.10 Vector distortions vs bit rate for compressed AVIRIS Cuprite Mine ’97 data. Radiometric distortion (a) VRMSE; spectral distortion (b) SAM
associated with low responses from the spectrometer. This explains why the bit rate plots of Fig. 6.9b never fall below 1 bit/pixel per band. Some of the radiometric distortion measures defined in Sect. 6.1 have been calculated on the distorted hyperspectral pixel vectors achieved by decompressing the bit streams generated by the near-lossless 3D-RLP encoder, both MAD- and PMAD-bounded. VRMSEs of the vector data, both average (6.18) and maximum (6.19), are plotted in Fig. 6.10a as a function of the bit rate from the encoder. The MAD-bounded encoder obviously minimizes both the average and maximum of VRMSE, that is the Euclidean norm of the pixel error vector. A further advantage is that average and maximum VRMSE are very close to each other for all bit rates. The PMAD-bounded encoder is somewhat poorer: average VRMSE is comparable with that of the former, but peak VRMSE is far larger, due to the highsignal components that are coarsely quantized in order to minimize PMAD. Trivial results, not reported in the plots, are that MAD of the data cube (6.22) is exactly equal to the desired value, whereas the PMAD, being unconstrained, is higher. Symmetric results have been found by measuring PMAD on MAD-bounded and PMAD-bounded decoded data. As far as radiometric distortion is concerned, results are not surprising. Radiometric distortions measured on vectors are straightforwardly derived from those measured on scalar pixel values. The introduction of such spectral measurements as SAM (6.24) and SID (6.28) may overcome the rationale of distortion, as established in the signal/image processing community. Figure 6.10b shows spectral distortions between original and decompressed hyperspectral pixel vectors. The PMADbounded algorithm yields plots that lie in the middle between the corresponding ones produced by the MAD-bounded algorithm and are very close to each other too. Since the maximum SAM is a better clue of spectral quality of the decoded data than the average SAM may be, a likely conclusion would be that PMAD-bounded compression optimizes the spectral quality of the data, while MAD-bounded is superior in terms of radiometric quality. Furthermore, the maximum SAM
144
B. Aiazzi et al.
introduced by the P-MAD bounded logarithmic quantizer is lower than 0.2 ∘ for an average rate of 1 bit/pixel per vector component, i.e., CR¼16.
8.3
Discrimination of Materials
Analysis procedures of hyperspectral vectors are usually performed on reflectance data, especially when the goal is identification of materials by comparing their remotely sensed spectra with sample spectra extracted from reference spectral libraries. Whenever measured spectra are to be compared to laboratory spectra, the radiance data are converted into reflectances, e.g., by means of the following simplified formula: rðlÞ ¼
RðlÞ p IðlÞ TðlÞ
(6.36)
in which r(l) is reflectance, I(l) is solar irradiance on ground, T(l) is atmospheric transmittance, and R(l) is at-sensor radiance, all functions of wavelength l. Distortions introduced by compression on radiance will be amplified or attenuated depending on the values of the product I(l) T(l). So, spectral distortion, e.g., SAM, must be measured between reflectance pixel spectra, rather than radiance pixel spectra. Extensive results of spectral discrimination from compressed hyperspectral data have demonstrated that a SAM distortion lower than 0.5 ∘ has negligible impact on the capability of automated classifiers of identifying spectral signatures of materials [26]. As a matter of fact, uniform distortion allocation over the entire spectrum yields minimum angle error between original and decompressed vectors. Band-variable distortion allocation following the virtually lossless protocol does not minimize error. However, when decompressed radiance data are converted into reflectance data by means of (6.36), the distribution of errors in decompressed reflectance spectra becomes approximately flat with the wavelength. In fact, quantization step sizes for virtually lossless compression follow the trend in Fig. 6.5, while the product I(l) T (l) by which compression errors on radiance are multiplied to yield compression errors on reflectance follows an opposite trend. Eventually, virtually lossless compression is preferable to near-lossless compression at the same rate, at least for detection of materials using spectral libraries, because SAM of reflectances produced from virtually lossless data are in average lower than those coming from near lossless data.
6 Quality Issues for Compression of Hyperspectral Imagery. . .
145
9 Conclusions This chapter has pointed out that quality issues are crucial for compression of radiance data produced by imaging spectrometers. Most widely used lossless compression techniques may be possibly replaced by lossy, yet error-bounded, techniques (near-lossless compression). Both lossless and near-lossless compression techniques can be implemented as adaptive DPCM encoders, with minor differences between the former and the latter. The rationale that compression-induced distortion is more tolerable, i.e., less harmful, in those bands, in which the noise is higher, and vice-versa, constitutes the virtually lossless paradigm, which provides operational criteria to design quantization in DPCM schemes. The quantization step size of each band should be a fraction of the measured standard deviation of instrumental noise. For typical VNIR+SWIR (400–2500 nm) spectra, virtually lossless compression exploits a band-varying quantization, with step sizes approximately decaying with increasing wavelengths. Once decompressed radiance spectra are converted to reflectance by removing the contribution of solar irradiance and atmospheric transmittance, the distribution of compression-induced errors with the wavelength is roughly equalized. Hence, angular errors introduced by compression errors on radiances will be lower for virtually lossless compression than for near-lossless compression. This feature is expected to be useful for discrimination of materials from compressed data and spectral libraries.
References 1. Abrardo, A., Alparone, L., Bartolini, F.: Encoding-interleaved hierarchical interpolation for lossless image compression. Signal Processing 56(2), 321–328 (1997) 2. Aiazzi, B., Alba, P., Alparone, L., Baronti, S.: Lossless compression of multi/hyper-spectral imagery based on a 3-D fuzzy prediction. IEEE Trans. Geosci. Remote Sensing 37(5), 2287–2294 (1999) 3. Aiazzi, B., Alparone, L., Barducci, A., Baronti, S., Marcoionni, P., Pippi, I., Selva, M.: Noise modelling and estimation of hyperspectral data from airborne imaging spectrometers. Annals of Geophysics 41(1), 1–9 (2006) 4. Aiazzi, B., Alparone, L., Barducci, A., Baronti, S., Pippi, I.: Estimating noise and information of multispectral imagery. J. Optical Engin. 41(3), 656–668 (2002) 5. Aiazzi, B., Alparone, L., Baronti, S.: A reduced Laplacian pyramid for lossless and progressive image communication. IEEE Trans. Commun. 44(1), 18–22 (1996) 6. Aiazzi, B., Alparone, L., Baronti, S.: Near-lossless compression of 3-D optical data. IEEE Trans. Geosci. Remote Sensing 39(11), 2547–2557 (2001) 7. Aiazzi, B., Alparone, L., Baronti, S.: Context modeling for near-lossless image coding. IEEE Signal Processing Lett. 9(3), 77–80 (2002) 8. Aiazzi, B., Alparone, L., Baronti, S.: Fuzzy logic-based matching pursuits for lossless predictive coding of still images. IEEE Trans. Fuzzy Systems 10(4), 473–483 (2002) 9. Aiazzi, B., Alparone, L., Baronti, S.: Near-lossless image compression by relaxation-labelled prediction. Signal Processing 82(11), 1619–1631 (2002) 10. Aiazzi, B., Alparone, L., Baronti, S.: Lossless compression of hyperspectral images using multiband lookup tables. IEEE Signal Processing Lett. 16(6), 481–484 (2009)
146
B. Aiazzi et al.
11. Aiazzi, B., Alparone, L., Baronti, S., Lastri, C.: Crisp and fuzzy adaptive spectral predictions for lossless and near-lossless compression of hyperspectral imagery. IEEE Geosci. Remote Sens. Lett. 4(4), 532–536 (2007) 12. Aiazzi, B., Alparone, L., Baronti, S., Lotti, F.: Lossless image compression by quantization feedback in a content-driven enhanced Laplacian pyramid. IEEE Trans. Image Processing 6 (6), 831–843 (1997) 13. Aiazzi, B., Alparone, L., Baronti, S., Santurri, L.: Near-lossless compression of multi/ hyperspectral images based on a fuzzy-matching-pursuits interband prediction. In: S.B. Serpico (ed.) Image and Signal Processing for Remote Sensing VII, vol. 4541, pp. 252–263 (2002) 14. Alecu, A., Munteanu, A., Cornelis, J., Dewitte, S., Schelkens, P.: On the optimality of embedded deadzone scalar-quantizers for wavelet-based L-infinite-constrained image coding. IEEE Signal Processing Lett. 11(3), 367–370 (2004) 15. Alecu, A., Munteanu, A., Cornelis, J., Dewitte, S., Schelkens, P.: Wavelet-based scalable L-infinity-oriented compression. IEEE Trans Image Processing 15(9), 2499–2512 (2006) 16. Baraldi, A., Blonda, P.: A survey of fuzzy clustering algorithms for pattern recognition–Parts I and II. IEEE Trans. Syst. Man Cybern.–B 29(6), 778–800 (1999) 17. Benazza-Benyahia, A., Pesquet, J.C., Hamdi, M.: Vector-lifting schemes for lossless coding and progressive archival of multispectral images. IEEE Trans. Geosci. Remote Sensing 40(9), 2011–2024 (2002) 18. Bezdek, J.C.: Pattern Recognition with Fuzzy Objective Function Algorithm. Plenum Press, New York (1981) 19. Carpentieri, B., Weinberger, M.J., Seroussi, G.: Lossless compression of continuous-tone images. Proc. of the IEEE 88(11), 1797–1809 (2000) 20. Chang, C.I.: An information-theoretic approach to spectral variability, similarity, and discrimination for hyperspectral image analysis. IEEE Trans. Inform. Theory 46(5), 1927–1932 (2000) 21. Deng, G., Ye, H., Cahill, L.W.: Adaptive combination of linear predictors for lossless image compression. IEE Proc.-Sci. Meas. Technol. 147(6), 414–419 (2000) 22. Golchin, F., Paliwal, K.K.: Classified adaptive prediction and entropy coding for lossless coding of images. In: Proc. IEEE Int. Conf. on Image Processing, vol. III/III, pp. 110–113 (1997) 23. Huang, B., Sriraja, Y.: Lossless compression of hyperspectral imagery via lookup tables with predictor selection. In: L. Bruzzone (ed.) Proc. of SPIE, Image and Signal Processing for Remote Sensing XII, vol. 6365, pp. 63650L.1–63650L.8 (2006) 24. Jayant, N.S., Noll, P.: Digital Coding of Waveforms: Principles and Applications to Speech and Video. Prentice Hall, Englewood Cliffs, NJ (1984) 25. Ke, L., Marcellin, M.W.: Near-lossless image compression: minimum entropy, constrainederror DPCM. IEEE Trans. Image Processing 7(2), 225–228 (1998) 26. Keshava, N.: Distance metrics and band selection in hyperspectral processing with applications to material identification and spectral libraries. IEEE Trans. Geosci. Remote Sensing 42(7), 1552–1565 (2004) 27. Kiely, A.B., Klimesh, M.A.: Exploiting calibration-induced artifacts in lossless compression of hyperspectral imagery. IEEE Trans. Geosci. Remote Sensing 47(8), 2672–2678 (2009) 28. Klimesh, M.: Low-complexity adaptive lossless compression of hyperspectral imagery. In: Satellite Data Compression, Communication and Archiving II, Proc. SPIE, vol. 6300 pp. 63000N.1–63000N.9 (2006) 29. Lastri, C., Aiazzi, B., Alparone, L., Baronti, S.: Virtually lossless compression of astrophysical images. EURASIP Journal on Applied Signal Processing 2005(15), 2521–2535 (2005) 30. Magli, E., Olmo, G., Quacchio, E.: Optimized onboard lossless and near-lossless compression of hyperspectral data using CALIC. IEEE Geosci. Remote Sensing Lett. 1(1), 21–25 (2004) 31. Matsuda, I., Mori, H., Itoh, S.: Lossless coding of still images using minimum-rate predictors. In: Proc. IEEE Int. Conf. on Image Processing, vol. I/III, pp. 132–135 (2000)
6 Quality Issues for Compression of Hyperspectral Imagery. . .
147
32. Mielikainen, J.: Lossless compression of hyperspectral images using lookup tables. IEEE Signal Proc. Lett. 13(3), 157–160 (2006) 33. Mielikainen, J., Toivanen, P.: Clustered DPCM for the lossless compression of hyperspectral images. IEEE Trans. Geosci. Remote Sensing 41(12), 2943–2946 (2003) 34. Mielikainen, J., Toivanen, P., Kaarna, A.: Linear prediction in lossless compression of hyperspectral images. J. Optical Engin. 42(4), 1013–1017 (2003) 35. Penna, B., Tillo, T., Magli, E., Olmo, G.: Progressive 3-D coding of hyperspectral images based on JPEG 2000. IEEE Geosci. Remote Sensing Lett. 3(1), 125–129 (2006) 36. Pennebaker, W.B., Mitchell, J.L.: JPEG: Still Image Compression Standard. Van Nostrand Reinhold, New York (1993) 37. Ramabadran, T.V., Chen, K.: The use of contextual information in the reversible compression of medical images. IEEE Trans. Medical Imaging 11(2), 185–195 (1992) 38. Rao, A.K., Bhargava, S.: Multispectral data compression using bidirectional interband prediction. IEEE Trans. Geosci. Remote Sensing 34(2), 385–397 (1996) 39. Rao, K.K., Hwang, J.J.: Techniques and Standards for Image, Video, and Audio Coding. Prentice Hall, Engl. Cliffs, NJ (1996) 40. Reichel, J., Menegaz, G., Nadenau, M.J., Kunt, M.: Integer wavelet transform for embedded lossy to lossless image compression. IEEE Trans. Image Processing 10(3), 383–392 (2001) 41. Rice, R.F., Plaunt, J.R.: Adaptive variable-length coding for efficient compression of spacecraft television data. IEEE Trans. Commun. Technol. COM-19(6), 889–897 (1971) 42. Rizzo, F., Carpentieri, B., Motta, G., Storer, J.A.: Low-complexity lossless compression of hyperspectral imagery via linear prediction. IEEE Signal Processing Lett. 12(2), 138–141 (2005) 43. Roger, R.E., Cavenor, M.C.: Lossless compression of AVIRIS images. IEEE Trans. Image Processing 5(5), 713–719 (1996) 44. Said, A., Pearlman, W.A.: An image multiresolution representation for lossless and lossy compression. IEEE Trans. Image Processing 5(9), 1303–1310 (1996) 45. Tate, S.R.: Band ordering in lossless compression of multispectral images. IEEE Trans. Comput. 46(4), 477–483 (1997) 46. Taubman, D.S., Marcellin, M.W.: JPEG2000: Image compression fundamentals, standards and practice. Kluwer Academic Publishers, Dordrecht, The Netherlands (2001) 47. Wang, J., Zhang, K., Tang, S.: Spectral and spatial decorrelation of Landsat-TM data for lossless compression. IEEE Trans. Geosci. Remote Sensing 33(5), 1277–1285 (1995) 48. Weinberger, M.J., Rissanen, J.J., Arps, R.B.: Applications of universal context modeling to lossless compression of gray-scale images. IEEE Trans. Image Processing 5(4), 575–586 (1996) 49. Weinberger, M.J., Seroussi, G., Sapiro, G.: The LOCO-I lossless image compression algorithm: principles and standardization into JPEG-LS. IEEE Trans. Image Processing 9(8), 1309–1324 (2000) 50. Witten, I.H., Neal, R.M., Cleary, J.G.: Arithmetic coding for data compression. Commun. ACM 30, 520–540 (1987) 51. Wu, X., Bao, P.: L1 constrained high-fidelity image compression via adaptive context modeling. IEEE Trans. Image Processing 9(4), 536–542 (2000) 52. Wu, X., Memon, N.: Context-based, adaptive, lossless image coding. IEEE Trans. Commun. 45(4), 437–444 (1997) 53. Wu, X., Memon, N.: Context-based lossless interband compression–Extending CALIC. IEEE Trans. Image Processing 9(6), 994–1001 (2000)
Chapter 7
Ultraspectral Sounder Data Compression by the Prediction-Based Lower Triangular Transform Shih-Chieh Wei and Bormin Huang
Abstract The Karhunen–Loeve transform (KLT) is the optimal unitary transform that yields the maximum coding gain. The prediction-based lower triangular transform (PLT) features the same decorrelation and coding gain properties as KLT but with lower complexity. Unlike KLT, PLT has the perfect reconstruction property which allows its direct use for lossless compression. In this paper, we apply PLT to carry out lossless compression of the ultraspectral sounder data. The experiment on the standard ultraspectral test dataset of ten AIRS digital count granules shows that the PLT compression scheme compares favorably with JPEG-LS, JPEG2000, LUT, SPIHT, and CCSDS IDC 5/3.
1 Introduction Contemporary and future ultraspectral infrared sounders such as AIRS [1], CrIS [2], IASI [3] and GIFTS [4] represent a significant technical advancements in for environmental and meteorological prediction and monitoring. Given the large 3D volume of data obtained from high spectral and spatial observations, the use of effective data compression techniques will be beneficial for data transfer and storage. When the ultraspectral sounder data is used to retrieve geophysical parameters like the vertical profiles of atmospheric temperature, moisture and trace gases, the retrieval process involves solving the radiative transfer equation which is an ill-posed inverse problem and sensitive to the noise and error in the
S.-C. Wei (*) Department of Information Management, Tamkang University, Tamsui, Taiwan e-mail: [emailprotected] B. Huang Space Science and Engineering Center, University of Wisconsin, Madison, USA e-mail: [emailprotected] B. Huang (ed.), Satellite Data Compression, DOI 10.1007/978-1-4614-1183-3_7, # Springer Science+Business Media, LLC 2011
149
150
S.-C. Wei and B. Huang
data [5]. Thus, lossless or near lossless compression of the data is desired in order to avoid substantial degradation during retrieval. Past studies on lossless compression of the ultraspectral sounder data can be categorized into clustering-based, prediction-based and transform-based methods [6]. Since a higher correlation exists in spectral dimension than in spatial dimension of ultraspectral sounder data [6], there were studies on band reordering as preprocessing before compression [7, 8]. The Karhunen–Loeve transform (KLT), a.k.a. the principal component analysis (PCA) or the Hoteling transform, is the optimal unitary transform that yields the maximum coding gain. However, considering its data-dependent computational cost involving the eigenvectors from the input covariance matrix, KLT is often only used as a benchmark for performance comparison. Phoong and Lin [9] developed the prediction-based lower triangular transform (PLT) that features the same decorrelation and coding gain properties as KLT but with a lower design and implementational cost. They showed promising results when PLT was applied to lossy compression of 2D imagery and the AR(1) process [9]. However, the original PLT by Phoong et al. requires the input vector to be a blocked version of a scalar wide sense stationary (WSS) process. Weng et al. [10] proposed a generalized triangular decomposition (GTD) which allows the input vector to be a vector WSS process. GTD has the same coding gain as KLT and it includes KLT and PLT as special cases [10]. Furthermore, unlike KLT, PLT has a perfect reconstruction (PR) property which makes it usefull for lossless compression. Since PLT provides the same coding gain as KLT, with lower complexity and PR property, we were motivated to apply PLT to carry out lossless compression of the ultraspectral sounder data. The compression method consists of using PLT for spectral prediction, followed by the arithmetic coding. The PLT compression ratio will be compared with prediction-based methods like LUT [11] and JPEG-LS [12], and wavelet-transform-based methods like JPEG2000 [13], SPIHT [14] and CCSDS IDC 5/3 [15]. The rest of the paper is arranged as follows. Section 2 describes the ultraspectral sounder data used in this study. Section 3 introduces our compression scheme. Section 4 shows the compression results on the ultraspectral sounder data. Section 5 gives the conclusions.
2 Data The ultraspectral sounder data can be generated from either a Michelson interferometer (e.g. CrIS [2], IASI [3], GIFTS [4]) or a grating spectrometer (e.g. AIRS [1]). A standard ultraspectral sounder data set for compression is publicly available via anonymous ftp at ftp://ftp.ssec.wisc.edu/pub/bormin/Count/. It consists of ten granules, five daytime and five nighttime, selected from representative geographical regions of the Earth. Their locations, UTC times and local time adjustments are listed in Table 7.1. This standard ultraspectral sounder data set adopts the NASA EOS AIRS digital counts made on March 2, 2004. The AIRS data
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
151
includes 2,378 infrared channels in the 3.74–15.4 mm region of the spectrum. A day’s worth of AIRS data is divided into 240 granules, each of 6 min durations. Each granule consists of 135 scan lines containing 90 cross-track footprints per scan line; thus there are a total of 135 90 ¼ 12,150 footprints per granule. More information regarding the AIRS instrument may be acquired from the NASA AIRS web site at http://www-airs.jpl.nasa.gov. The digital count data ranges from 12 to 14 bits for different channels. Each channel is saved using its own bit depth. To make the selected data more generic to other ultraspectral sounders, 271 AIR-specific bad channels identified in the supplied AIRS infrared channel properties file are excluded. Each resulting granule is saved as a binary file, arranged as 2,107 channels, 135 scan lines, and 90 pixels for each scan line. Figure 7.1 shows the AIRS digital counts at wavenumber 800.01 cm 1 for some of the ten selected granules. In these granules, coast lines are depicted by solid curves and multiple clouds at various altitudes are shown as different shades of pixels. Table 7.1 Ten selected airs granules for study of ultraspectral sounder data compression Granule number UTC time Local time adjustment Location Granule 9 00:53:31 UTC 12 H Pacific Ocean, daytime Granule 16 01:35:31 UTC +2 H Europe, nighttime Granule 60 05:59:31 UTC +7 H Asia, daytime Granule 82 08:11:31 UTC 5H North America, nighttime Granule 120 11:59:31 UTC 10 H Antarctica, nighttime Granule 126 12:35:31 UTC 0H Africa, daytime Granule 129 12:53:31 UTC 2H Arctic, daytime Granule 151 15:05:31 UTC +11 H Australia, nighttime Granule 182 18:11:31 UTC +8 H Asia, nighttime Granule 193 19:17:31 UTC 7H North America, daytime
Fig. 7.1 The ten selected AIRS digital count granules at wavenumber 800.01 cm 2004
1
on March 2,
152
Fig. 7.1 (continued)
S.-C. Wei and B. Huang
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
153
Fig. 7.1 (continued)
3 The Compression Scheme For ultraspectral sounder data, the spectral correlation is generally much stronger than the spatial correlation [6]. To de-correlate the spectral dependency, a linear prediction with a fixed number of predictors has been used [16]. The prediction-based lower triangular transform (PLT) is also based on linear prediction but uses as many orders of predictors as possible. However, without doing linear regression on each order, PLT can be directly computed by the LDU matrix decomposition [9]. Figure 7.2 gives the schematic of the PLT transform coding scheme. Let x(n) be a sequence of scalar observation signals. In order to exploit the correlation in the sequence, M consecutive scalar signals are grouped together to form a blocked version of vector signal ~ xðtÞ. Each vector signal ~ xðtÞ then goes through a PLT transform T to obtain a vector of transform coefficients ~ yðtÞ. The transform coefficients are often quantized and have smaller variances than the original signal for storage or transfer purposes. Finally, when necessary, an inverse PLT transform T 1 is applied to the transform coefficients to restore the original signal. To compute the PLT transform and the inverse transform, the statistic of the source signal are required. Suppose a total of N M samples are collected, the N M signal matrix X can be expressed as: 0 1 xð1Þ xðM þ 1Þ xð2M þ 1Þ xððN 1ÞM þ 1Þ B xð2Þ xðM þ 2Þ xð2M þ 2Þ xððN 1ÞM þ 2Þ C B C B C B xð3Þ xðM þ 3Þ xð2M þ 3Þ xððN 1ÞM þ 3Þ C X¼B C C B . .. .. .. .. C B . . . . . A @ . xðMÞ xð2MÞ xð3MÞ xðNMÞ ¼ ð~ x(1) ~ x(2) ~ x(3) ~ x(N))
154
S.-C. Wei and B. Huang
Fig. 7.2 The schematic of the original PLT transform coding scheme
x(1)
y(1)
x(2)
y(2)
x(3)
T
x(M)
y(3)
x(1) x(2) -1
T
y(M)
x(3)
x(M)
Assume that Rx is the autocorrelation matrix of the observation signal x(n) with order M, by the LDU decomposition of Rx, the PLT transform for the signal X can be computed as follows: P ¼ L 1 with Rx ¼ LDU where L is the lower triangular matrix with all diagonal entries equal to 1, D is a diagonal matrix, U is the upper triangular matrix with all diagonal entries equal to 1, and P is the desired prediction-based lower triangular (PLT) transform T which can be computed by the matrix inversion of L. Since Rx is symmetric and the LDU matrix decomposition is unique [17], we have LDU ¼ Rx ¼ RxT ¼ U T D LT or U ¼ LT. Let y(n) be the coefficient of the transform P on x(n). y(n) will be the error between the original signal x(n) and the prediction based on previous M 1 signals, i.e. x(n 1) through x(n M + 1). Let Ry be the autocorrelation matrix of the prediction error y(n). By the transform property Ry ¼ P Rx(M) PT and the LDU decomposition of Rx ¼ LDU, Ry can be shown to be diagonal as follows [9]: Ry ¼ PRx PT ¼ PðLDUÞPT ¼ PP 1 DðPT Þ 1 PT ¼ D It means that the prediction error is de-correlated by the initial selection of P ¼ L 1. The entries on the diagonal of D will be the variance of the prediction error of orders 0 through M 1. Since the lower triangular transform P consists of all linear predictors below order M, computation of P can be carried out by the Levinson-Durbin algorithm [18] which has lower complexity than the LDU decomposition when the scalar signal satisfies the wide sense stationary (WSS) property. For compression of the ultraspectral sounder data, each granule is considered to consist of N ¼ ns observations at different locations with each observation containing M ¼ nc channels of data. Specifically, we select the prediction order M as the number of channels nc to fully exploit the spectral correlation in reducing the prediction error. That is, all previous available channels are used to predict the current channel. Figure 7.3 gives the schematic of the PLT transform coding model for our compression scheme. Let X ¼ [x1 x2 x3 . . . xnc]T be the original mean-subtracted ultraspectral sounder data consisting of nc channels by ns observations, Y ¼ [y1 y2 y3 . . . ync]T be the nc x
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . . Fig. 7.3 The schematic of the PLT transform coding model for our compression scheme
155
x1
y1
x1
x2
y2
x2
x3
T-1
y3
T
xnc
x3
xnc
ync
ns transform coefficients or prediction errors, and P be the prediction-based lower triangular transform. Here X, Y, and P can be written as follows: 1 0 1 0 x1 ð1Þ x1 ð2Þ x1 ð3Þ x1 ðns Þ ~ x1 B C B x ð1Þ x ð2Þ x ð3Þ x ðn Þ C x2 C B 2 B~ 2 2 2 s C C B C B C B x3 ð1Þ x3 ð2Þ x3 ð3Þ x3 ðns Þ C B~ x 3 X¼B C¼B C B . C B . .. C .. .. .. B . C B . C . A . . . @ . A @ . ~ xnc xnc ð1Þ xnc ð2Þ xnc ð3Þ xnc ðns Þ ¼ð~ x(1) ~ x(2) ~ x(3) ~ x(ns)) 0 1 0 ~ y1 y1 ð1Þ y1 ð2Þ C B~ B B y2 C B y2 ð1Þ y2 ð2Þ B C B B y3 C B y3 ð1Þ y3 ð2Þ Y ¼B~ C¼B B . C B . .. B . C B . . @ . A @ . ~ ync ync ð1Þ ync ð2Þ
y1 ð3Þ
y1 ðns Þ
y2 ð3Þ y3 ð3Þ .. . ync ð3Þ
.. .
y2 ðns Þ C C C y3 ðns Þ C C .. C C . A ync ðns Þ
¼ ð~ y(1) ~ y(2) ~ y(3) ~ y(ns)) 0 1 0 0 B p1;0 1 0 B B p2;1 1 P ¼ B p2;0 B .. .. .. @ . . . pnc 1;0 pnc 1;1 pnc 1;2
1 0 0C C 0C C .. C .. . .A 1
1
Then the transform coefficient or prediction error Y can be computed by Y ¼ PX or y1 ¼ x1 y2 ¼ x2 y3 ¼ x3
x^2 ¼ x2 þ p1;0 x1 x^3 ¼ x3 þ p2;0 x1 þ p2;1 x2
::: ync ¼ xnc
x^nc
¼ xnc þ pnc
1;0 x1
þ pn c
1;1 x2
þ pnc
1;2 x3
þ ::: þ pnc
1;nc 2 xnc 1
where x^m is the prediction of channel m by use of the linear combination of all previous m 1 channels. A similar result can be obtained for the inverse transform
156
S.-C. Wei and B. Huang
S ¼ P 1 which is derived from L in the LDU decomposition. Let the inverse transform S, which is also a lower triangular matrix, be of the form 1 0 1 0 0 0 B s1;0 1 0 0C C B B s2;0 s 1 0C 2;1 S¼B C B .. .. C .. .. .. @ . . .A . . snc 1;0 snc 1;1 snc 1;2 1 Then the original signal X can be computed by X ¼ SY or x1 ¼ y1 x2 ¼ y2 þ x^2 ¼ y2 þ s1;0 y1 x3 ¼ y3 þ x^3 ¼ y3 þ s2;0 y1 þ s2;1 y2 ::: xnc ¼ ync þ x^nc ¼ y n c þ sn c
1;0 y1
þ snc
1;1 y2
þ snc
1;2 y3
þ ::: þ snc
1;nc 2 ync 1
From the above formulation, it can be seen that either the transform P or the inverse transform S alone can be used to compress the signal X into the prediction error Y in order to reconstruct the original X from Y. However, to minimize the amount of data in transfer, quantization on transform kernels P or S and quantization on the prediction error Y are required. Furthermore, perfect reconstruction is required for our lossless compression application. To meet these requirements, a minimum noise ladder-based structure with the perfect reconstruction property is adopted [9]. When the transform P is used for both encoding and decoding, a minimum noise ladder-based structure for encoding follows. Note that the flooring function is chosen for quantization. Both X and Y will be integers. y1 ¼ x1 y2 ¼ x2 þ floorðp1;0 x1 Þ y3 ¼ x3 þ floorðp2;0 x1 þ p2;1 x2 Þ ::: ync ¼ xnc þ floorðpnc
1;0 x1
þ pnc
1;1 x2
þ pn c
1;2 x3
þ ::: þ pnc
1;nc 2 xnc 1 Þ
A corresponding minimum noise ladder-based structure using P for decoding follows: x1 ¼ y1 x2 ¼ y2
floorðp1;0 x1 Þ
x3 ¼ y3 :::
floorðp2;0 x1 þ p2;1 x2 Þ
xnc ¼ ync
floorðpnc
1;0 x1
þ pnc
1;1 x2
þ pn c
1;2 x3
þ ::: þ pnc
1;nc 2 xnc 1 Þ:
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
x1 x2 x3
pnc-1,0
p2,0 p1,0 floor
p2,1
pnc-1,1
floor
pnc-1,2
y1 y2
x1
-p1,0 -p2,0 floor
-pnc-1,0
-p2,1
x2
-pnc-1,1
x3
floor
y3
-pnc-1,2
pnc-1,nc-2 xnc
157
-pnc-1,nc-2 ync
floor
floor
xnc
Fig. 7.4 A minimum noise ladder-based structure of PLT based on the transform P and the flooring quantizer. It has the perfect reconstruction property for lossless compression
Similarly when the inverse transform S is used for both encoding and decoding, a minimum noise ladder-based structure with the perfect reconstruction property follows. Note that the ceiling function is chosen here for quantization. y1 ¼ x1 y2 ¼ x2
ceilðs1;0 y1 Þ
y3 ¼ x3 :::
ceilðs2;0 y1 þ s2;1 y2 Þ
ync ¼ xnc
ceilðsnc
1;0 y1
þ snc
1;1 y2
þ sn c
1;2 y3
þ ::: þ snc
1;nc 2 ync 1 Þ
A corresponding minimum noise ladder-based structure using S for decoding follows. x1 ¼ y1 x2 ¼ y2 þ ceilðs1;0 y1 Þ x3 ¼ y3 þ ceilðs2;0 y1 þ s2;1 y2 Þ ::: xnc ¼ ync þ ceilðsnc
1;0 y1
þ snc
1;1 y2
þ sn c
1;2 y3
þ ::: þ snc
1;nc 2 ync 1 Þ
Figure 7.4 shows the diagram for a minimum noise ladder-based structure of PLT for encoding and decoding where the transform P and the flooring quantizer are used. Similarly, Fig. 7.5 shows the diagram for a minimum noise ladder-based structure of PLT for encoding and decoding where the inverse transform S and the ceiling quantizer are used. Note that in both designs, all inputs to the multipliers of the transform kernel are quantized values so that these same values can be used in reconstruction to restore the original values. To save storage space, quantization is not only applied to
158
x1 x2
S.-C. Wei and B. Huang
-s1,0 -s2,0 ceil
x3
-s2,1
y1
-snc-1,0
y2
-snc-1,1
y3
ceil
-snc-1,2
ceil
snc-1,1
ceil
s2,1 ceil
snc-1,2
x2 x3
snc-1,nc-2
-snc-1,nc-2 xnc
x1
s2,0 s1,0
snc-1,0
ync
xnc
ceil
Fig. 7.5 A minimum noise ladder-based structure of PLT based on the inverse transform S and the ceiling quantizer. It also has the perfect reconstruction property for lossless compression
Original Granule
Mean Subtraction
Mean
Mean Subtracted Granule
Compute Transform Kernel
Transform Kernel
Quantizer
PLT Transform
Quantized Transform Kernel
Prediction Error
Arithmetic Coding
Compressed Granule
Fig. 7.6 The data flow diagram of the PLT compression scheme
the transform coefficients but also to the transform kernel. Furthermore, the traditional standard scalar or vector quantization techniques, which require transmission of the codebook, are not used. Instead, a simple ceiling or flooring on a specified number of decimal places is used for the codebook free quantization. With the minimum noise ladder-based structure being used, only the quantized transform kernel and the quantized prediction error need to be sent for restoration. Both then go through the arithmetic coder [19] in order to enhance the compression ratio. The data flow diagrams of our PLT compression and decompression schemes are shown in Figs. 7.6 and 7.7 respectively. As the original granule in Fig. 7.6 is meansubtracted, in addition to the compressed transform kernel and the compressed prediction error, a mean vector of nc 1 has to be sent to the decoder in Fig. 7.7 to recover the original granule.
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
Compressed Granule
Arithmetic Decoding
Quantized Transform Kernel
PLT Inverse Transform
159
+
Recovered Granule
Prediction Error Mean
Fig. 7.7 The data flow diagram of the PLT decompression scheme
4 Results The standard test data of ten AIRS ultraspectral digital count granules are tested. As performance index, the compression ratio and the bit rate are used. Both minimum noise ladder-based structures of PLT based on the transform P (PLT-P) and the inverse transform S (PLT-S) are tested. The quantizer used in PLT-S is tuned to keep the precision of the multipliers in transform S to two decimal places while the quantizer used in PLT-P is tuned to keep the precision of the multipliers in transform P to three decimal places. As the transformation process uses ceiling or flooring quantizers, the transform coefficient or the prediction error is an integer. A context-based arithmetic coder [20] is used to encode the transform kernel and the prediction error. The ten plots in Fig. 7.8 show the variances of channels for the ten AIRS test granules before and after the PLT transform. The dotted curve for X is the variances of channels before the PLT transform. The solid curve for Y is the variances of channels after the PLT transform. Throughout the ten semilog plots for the ten granules, it can be seen that the variances of Y are significantly lower than those of X at most channels. The lower variances of the transform coefficients Y promise less bits in need for compression. Figure 7.9 shows the amount of energy compaction in terms of the AM/GM ratio for the ten granules. The AM of a signal is the arithmetic mean of the variances of its subbands. The GM of a signal is the geometric mean of the variances of its subbands. It is shown that the ratio of AM to GM will always be greater than or equal to 1 [21]. AM will be equal to GM only when the variances of all subbands are equal. When a signal has a higher AM/GM ratio, there is a higher amount of energy compaction in its subband signals and the signal is good for data compression [22]. In Fig. 7.9, the AM/GM ratio is consistently much higher for the transformed coefficient Y than for the original signal X. This means that the energy compaction is much higher after the PLT transform than before the transform. Figure 7.10 shows the coding gain of the PLT transform. The coding gain of a transform coding T is defined as the mean squared reconstruction error in PCM coding divided by the mean squared reconstruction error in the transform coding T [23]. The reconstruction error is the absolute difference between the reconstructed signal and the original signal. In fact, the coding gain of the PLT transform can be
160
S.-C. Wei and B. Huang
Fig. 7.8 The variances of channels X and Y before and after the PLT transform respectively for the ten AIRS granules in the dataset
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
Fig. 7.8 (continued)
161
162
Fig. 7.8 (continued)
S.-C. Wei and B. Huang
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
Fig. 7.8 (continued)
163
164
Fig. 7.8 (continued)
S.-C. Wei and B. Huang
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
165
14 12
AM/GM
10 8 X Y
6 4 2 0 9
16
60
82
120 126 129 151 182 193 Granule No.
Fig. 7.9 The amount of energy compaction of the original signal X and the transformed signal Y as measured by the ratio of AM/GM for the ten AIRS granules in the dataset
100000
Coding Gain
10000 1000 100 10 1
9
16
60
82
120 126 129
151 182 193
Granule No.
Fig. 7.10 The coding gain of the PLT transform for the ten AIRS granules in the dataset. The PLT transform has the same maximum coding gain as the KLT transform
computed by dividing the arithmetic mean of the subband variances (AM) of the original signal by the geometric mean of the subband variances (GM) of the transform coefficient [9]. In Fig. 7.10, the PLT transform has achieved the maximum coding gain, same as the KLT transform. Figure 7.11 shows the compression ratio of the ten tested granules. The compression ratio is defined as the size of the original file divided by the size of the compressed file. For the AIRS granules, an original file containing 135902107 samples takes up about 41.2 MB. A higher compression ratio denotes a better compression result. The figure shows that PLT-S is slightly better than PLT-P.
166
S.-C. Wei and B. Huang
PLT-S
PLT-P
JPEG-LS
JPEG2000
LUT
SPIHT
CCSDS IDC 5/3
Compression Ratio
3.0
2.5
2.0
1.5
1.0
9
16
60
82
120
126
129
151
182
193
Granule Number
Fig. 7.11 The compression ratios of our compression scheme using the transform P (PLT-P) and the inverse transform S (PLT-S) in comparison with JPEG-LS, JPEG2000, LUT, SPIHT, and CCSDS IDC 5/3 on the ten ultraspectral granules
Table 7.2 Bit rates for our scheme using the transform P (PLT-P) and the inverse transform S (PLT-S) compared with JPEG-LS, JPEG2000, LUT, SPIHT, and CCSDS IDC 5/3 on the ten ultraspectral granules Granule number PLT-S PLT-P JPEG-LS JPEG2000 LUT SPIHT CCSDS IDC 5/3 9 4.21 4.26 5.24 5.46 5.63 6.48 6.45 16 4.23 4.19 5.14 5.26 5.39 6.29 6.29 60 4.32 4.34 5.37 5.75 5.73 6.86 6.82 82 4.22 4.19 5.00 5.06 5.30 6.08 6.05 120 4.28 4.27 5.20 5.68 5.65 6.71 6.68 126 4.31 4.37 5.37 5.60 5.68 6.68 6.68 129 4.15 4.17 5.00 5.14 6.11 6.17 6.14 151 4.42 4.41 5.28 5.81 5.81 6.82 6.82 182 4.44 4.47 5.44 5.65 5.94 6.68 6.68 193 4.30 4.33 5.35 5.81 5.81 6.86 6.86
Moreover both PLT-S and PLT-P compare favorably with the prediction-based methods like LUT [11], JPEG-LS [12], the wavelet-transform-based methods like JPEG2000 [13], SPIHT [14] and the CCSDS recommendation for Image Data Compression (IDC)-5/3 scheme [15, 24] in all 10 granules [24]. The alternative compression result in terms of the bit rate is summarized in Table 7.2. The bit rate is defined as the number of bits used per sample. A lower bit rate means that fewer bits are required to encode a sample and therefore denotes
7 Ultraspectral Sounder Data Compression by the Prediction-Based. . .
167
better result of compression. For the AIRS granules, a sample takes up 12–14 bits depending on its channel. In average, the bit rate of the original file is 12.9 bits per samples.
5 Conclusions The ultraspectral sounder data is characterized by its huge size and low tolerance of noise and error. Use of lossless or near lossless compression is thus desired for data transfer and storage. There have been prior works on lossless compression of ultraspectral sounder data which can be categorized into clustering-based, prediction-based, and transformation-based methods. In this work, a transformationbased method using the prediction-based lower triangular transform (PLT) is proposed. In our formulation, the PLT can use the P transform (PLT-P) or the S transform (PLT-S) for compression. To save space, a simple codebook-free quantization is applied to the transform kernel and the transform coefficient which is the prediction error for the PLT transform. Due to the minimum noise ladder structure in the PLT design, we can fully recover the granule from the quantized transform kernel and coefficients. To enhance the compression ratio, both the quantized transform kernel and the quantized transform coefficients are fed to the arithmetic coder for entropy coding. In terms of the compression ratio, the result shows that the PLT compression scheme compares favorably with the prediction-based methods like LUT and JPEG-LS, the wavelet-transform-based methods like JPEG2000, SPIHT and CCSDS IDC 5/3 on all the ten standard AIRS granules. However, one disadvantage of the method is that intensive CPU computation is required by our PLT compression scheme. Therefore, future work includes reduction of the CPU time requirements by adopting parallel computation.
References 1. H. H. Aumann and L. Strow, “AIRS, the first hyper-spectral infrared sounder for operational weather forecasting,” Proc. of IEEE Aerosp. Conf., 4, pp. 1683–1692, 2001. 2. H. J. Bloom, “The Cross-track Infrared Sounder (CrIS): a sensor for operational meteorological remote sensing,” Proc. of the 2001 Int. Geosci. and Remote Sens. Symp., pp. 1341–1343, 2001. 3. T. Phulpin, F. Cayla, G. Chalon, D. Diebel, and D. Schlussel, “IASI onboard Metop: Project status and scientific preparation,” 12th Int. TOVS Study Conf., Lorne, Victoria, Australia, pp. 234–243, 2002. 4. W. L. Smith, F. W. Harrison, D. E. Hinton, H. E. Revercomb, G. E. Bingham, R. Petersen, and J. C. Dodge, “GIFTS – the precursor geostationary satellite component of the future Earth Observing System,” Proc. of the 2002 Int. Geosci. and Remote Sens. Symp., 1, pp. 357–361, 2002. 5. B. Huang, W. L. Smith, H.-L. Huang, and H. M. Woolf, “Comparison of linear forms of the radiative transfer equation with analytic Jacobians”, Appl. Optics, vol. 41, no. 21, pp. 4209–4219, 2002.
168
S.-C. Wei and B. Huang
6. B. Huang, A. Ahuja, and H.-L. Huang, “Lossless compression of ultraspectral sounder data,” Hyperspectral Data Compression, G. Motta and J. Storer; Eds., Springer-Verlag, pp. 75–106, 2005. 7. P. Toivanen, O. Kubasova, and J. Mielikainen, “Correlation-based band-ordering heuristic for lossless compression of hyperspectral sounder data”, IEEE Geosci. Remote Sens. Lett., vol. 2, no. 1, pp.50–54, 2005. 8. B. Huang, A. Ahuja, H.-L. Huang, T. J. Schmit, and R. W. Heymann, “Lossless compression of 3D hyperspectral sounding data using context-based adaptive lossless image codec with biasadjusted reordering,” Optical Engineering, vol. 43, no. 9, pp. 2071–2079, 2004. 9. S.-M. Phoong and Y.-P. Lin, “Prediction-based lower triangular transform,” IEEE Trans. Signal Processing, vol. 48, no. 7, pp. 1947–1955, 2000. 10. C.-C. Weng, C.-Y. Chen and P. P. Vaidyanathan, “Generalized triangular decomposition in transform coding,” IEEE Trans. Signal Processing, vol. 58, no. 2, pp. 566–574, 2010. 11. B. Huang, and Y. Sriraja, “Lossless compression of hyperspectral imagery via lookup tables with predictor selection,” Proc. SPIE, vol. 6365, pp.63650L.1, 2006. 12. ISO/IEC 14495–1 and ITU Recommendation T.87, “Information Technology – lossless and near-lossless compression of continuous-tone still images,” 1999. 13. D. S. Taubman and M. W. Marcellin, JPEG2000: Image compression fundamentals, standards, and practice, 2002. 14. A. Said, and W. A. Pearlman, “A new, fast, and efficient image codec based on set partitioning in hierarchical trees,” IEEE Trans. Circuits. Sys. Video Tech., vol. 6, pp. 243–250, June 1996. 15. CCSDS, “Consultative Committee for Space Data Systems,” http://www.ccsds.org. 16. B. Huang, A. Ahuja, H.-L. Huang, T.J. Schmit, R.W. Heymann, “Fast precomputed VQ with optimal bit allocation for lossless compression of ultraspectral sounder data”, Proc. IEEE Data Comp. Conf., pp. 408–417, March 2005. 17. G. H. Golub and C. F. V. Loan, Matrix computations, John Hopkins University Press, 1996. 18. A. Gersho and R. M. Gray, Vector quantization and signal compression, Kluwer Academic Publishers, 1992. 19. I. Witten, R. Neal, and J. Cleary, “Arithmetic coding for data compression”, Comm. ACM, vol. 30, no. 6, pp. 520–540, June 1987. 20. M. R. Nelson, “Arithmetic coding and statistical modeling”, Dr. Dobb’s Journal, pp. 16–29, February 1991. 21. Y. You, Audio coding- theories and applications, Springer, 2010. 22. K. Sayood, Introduction to data compression, 2nd Ed., Morgan Kaufmann Publishers, 2000. 23. N. S. Jayant and P. Noll, Digital coding of waveforms- principles and applications to speech and video, Prentice Hall, 1984. 24. J. Serra-Sagrista, F. Garcia, J. Minguillon, D. Megias, B. Huang, and A. Ahuja, “Wavelet lossless compression of ultraspectral sounder data,” Proc. Int. Geosci. Rem. Sens. Symp., vol. 1, pp. 148–151, July 2005.
Chapter 8
Lookup-Table Based Hyperspectral Data Compression Jarno Mielikainen
Abstract This chapter gives an overview of the lookup table (LUT) based lossless compression methods for hyperspectral images. The LUT method searches the previous band for a pixel of equal value to the pixel co-located to the one to be coded. The pixel in the same position as the obtained pixel in the current band is used as the predictor. Lookup tables are used to speed up the search. Variants of the LUT method include predictor guided LUT method and multiband lookup tables.
1 Introduction Hyperspectral imagers produce enormous data volumes. Thus, a lot of effort has been spent to research more efficient ways to compress hyperspectral images. Three different types of compression modalities for hyperspectral images can be defined. Lossy compression achieves the lowest bit rate among the three modalities. It does not bind the difference between each reconstructed pixel and the original pixel. Instead, the reconstructed image is required to be similar to the original image on mean-squared error sense. Near lossless compression bounds the absolute difference between each reconstructed pixel and the original pixel by a predefined constant. Lossless compression requires the exact original image to be reconstructed from the compressed data. Since lossless compression techniques involve no loss of information they are used for applications that cannot tolerate any difference between the original and reconstructed data. In hyperspectral images the interband correlation is much stronger than the intraband correlation. Thus, interband correlation must be utilized for maximal compression performance. Transform-based and vector-quantization-based methods have not been able to achieve state-of-the-art lossless compression results
J. Mielikainen (*) School of Electrical and Electronic Engineering, Yonsei University, Seoul, South Korea e-mail: [emailprotected] B. Huang (ed.), Satellite Data Compression, DOI 10.1007/978-1-4614-1183-3_8, # Springer Science+Business Media, LLC 2011
169
170
J. Mielikainen
for hyperspectral images. Therefore, lossless compression of hyperspectral data is performed by using prediction-based approaches. However, there have been some studies on transform-based [1–3] and vector-quantization based [4–6] methods. Vector quantization is an asymmetric compression method; compression is much more computationally intensive than decompression. On the other hand, transformbased methods have been more successful in lossy compression than lossless compression. Prediction based methods for lossless compression of hyperspectral images can be seen as consisting of three steps: 1. Band ordering. 2. Modeling extracting information on the redundancy of the data and describing this redundancy in the form of a model. 3. Coding describes the model and how it differs from the data using a binary alphabet. The problem of optimal band ordering for hyperspectral image compression has been solved in [7]. Optimal band reordering is achieved by computing a minimum spanning tree for a directed graph containing the sizes of the encoded residual bands. A correlation-based heuristic for estimating the optimal order was proposed in [8]. Another prediction method based on reordering was introduced in [9]. However, in this chapter, all the experiments are performed using natural ordering of the bands to facilitate comparisons to the other methods in the literature. In this chapter, we concentrate on lookup table (LUT) based approaches to modeling and we are will gives an overview of LUT based lossless compression methods for hyperspectral images. This chapter is organized as follows. In Sect. 2 we will present a short review of previous work in lossless compression of hyperspectral images. Section 3 presents basic LUT method. In Sect. 4 predictor guided LUT is described. Use of a quantized index in LUT method is discussed in Sect. 5. Multiband generalization of LUT method is presented in Sect. 6. Experiments results are shown in Sect. 7. Finally, conclusions are drawn in Sect. 8.
2 Lossless Compression of Hyperspectral Images Previous approaches to lossless compression of hyperspectral images include A1, which is one of three distributed source coding algorithms proposed in [10]. It focuses on coding efficiency and the other two algorithms proposed in [10] are more focused on error-resiliency. The A1 algorithm independently encodes nonoverlapped blocks of 16 16 samples in each band. This independency makes it easy to parallelize the algorithm. The first block of each band is transmitted uncompressed. The pixel values are predicted by a linear prediction that utilizes pixel value in previous bans, the average pixel values of both the current block and the co-located block in the previous band. Instead of sending prediction parameters
8 Lookup-Table Based Hyperspectral Data Compression
171
to decoder they are guessed by the decoder. For each guess the pixels of the block are reconstructed and the Cyclic Redundancy Check (CRC) is computed. Once CRC matches the one included in the compressed file, the process terminates. The FL algorithm [11] employs the previous band for prediction and adapts the predictor coefficients using recursive estimation. The BG block-based compression algorithm [12] employs a simple block-based predictor followed by an adaptive Golomb code. IP3 (third-order interband predictor) [13] method takes advantage of spatial data correlation and derives spectral domain predictor using Wiener filtering. They also employed a special backward pixel search (BPS) module for calibrated image data. Clustered differential pulse code modulation (C-DPCM) [14] method partitions spectral vectors into clusters and then applies a separate least-squares optimized linear predictor to each cluster of each band. The method can be seen as an extension of the vector quantization method in [5]. However, the quantization step of [5] is omitted. In [15], another approach using clustering was presented. The causal neighborhoods of each pixel are clustered using fuzzy-c-means clustering. For each of the clusters, an optimal linear predictor is computed from the values, the membership degrees of which exceed a threshold. The final estimate is computed as a weighted sum of the predictors, where the weights are the membership degrees. The Spectral Fuzzy Matching Pursuits (S-FMP) method exploits a purely spectral prediction. In the same paper, a method called Spectral RelaxationLabeled Prediction (S-RLP) was also proposed. The method partitions image bands into blocks, and a predictor, out of a set of predictors, is selected for prediction. A method based on Context-Adaptive Lossless Image Coding (CALIC), which is called 3-D CALIC [28], switches between intra- and interband prediction modes based on the strength of the correlation between the consecutive bands. In multiband CALIC (M-CALIC) method [16], the prediction estimate is performed using two pixels in the previous bands in the same spatial position as the current pixel. The prediction coefficients are computed using an offline procedure on training data. An adaptive least squares optimized prediction technique called Spectrumoriented Least SQuares (SLSQ) was presented in [17]. The prediction technique used is the same as the one in [18], but a more advanced entropy coder was used. The predictor is optimized for each pixel and each band in a causal neighborhood of the current pixel. SLSQ-HEU uses a heuristic to select between the intra- and interband compression modes. Also, an optimal method for inter-/intracoding mode selection called SLSQ-OPT was presented. Selecting between a Correlation-based Conditional Average Prediction (CCAP) and a lossless JPEG was proposed in [19]. The selection is based on a correlation coefficient for contexts. The CCAP estimate is a sample mean of pixels corresponding to the current pixel in contexts that match the current pixel context. BH [20] is a block-based compressor. Each band of the input image is divided into square blocks. Next, the blocks are predicted based on the corresponding block in the previous band. Nonlinear Prediction for Hyperspectral Images (NPHI) [21] predicts the pixel in the current band based on the information in the causal context in the current band and pixels colocated in the reference band. NPHI was also extended
172
J. Mielikainen
into an edge-based technique, called the Edge-based Prediction for Hyperspectral Images, which classifies the pixels into edge and nonedge pixels. Each pixel is then predicted using information from pixels in the same pixel class within the context. In [23], a method called KSP, which employs a Kalman filter in the prediction stage, was proposed.
3 LUT Method The LUT method [22] makes a prediction of the current pixel px,y,z (xth row, yth column, and zth band) using all the causal pixels in the current and previous band. LUT method is based on the idea of Nearest Neighbor (NN) search. The NN procedure searches for the nearest neighbor in the previous band that has the same pixel value as the pixel located in the same spatial position as the current pixel in the previous band px,y,z 1. The search is performed in reverse raster-scan order. First, a pixel value equal to px,y,z 1 is searched. If an equal valued pixel is found at position (x’,y’,z 1), then estimated pixel is predicted to have the same value as the pixel in the same position as obtained pixel in the current band px’,y’,z. Otherwise, the estimated pixel value is equal to the pixel value in the previous band px,y,z 1. LUT method accelerates NN method by replacing time consuming search procedure with a lookup table operation, which uses the pixel co-located in the previous band as an index in the lookup table. The lookup table returns the nearest matching pixel. An example illustrating the search process is shown in Figs. 8.1–8.3. The example uses two consecutive image bands, which have 38 pixels each. The previous band (band number 1) and current band (band number 2) are shown in Figs. 8.1 and 8.2, respectively. The corresponding lookup table is shown in Fig. 8.3. In the example, pixel p3,8,2 ¼ 325 is the current pixel to be predicted in the current band. The causal pixels in the previous band are searched to find a match for the colocated pixel p3,8,1 ¼ 315. Both current pixel and its co-located pixel have yellow background in Figs. 8.2 and 8.1, respectively. Three matches (green background) are returned. The pixel value in the current band that is present at the nearest matching location, p2,6,1 ¼ 315, is used as the predictor for p’3,8,2 ¼ p2,6,2 ¼ 332. A time-consuming search was avoided because the lookup table directly returned the predictor value.
4 Predictor Guided LUT Method In the LUT method the nearest matching pixel value might be not be as good of a match as many other matching pixels. In the previous example the pixels in the current band corresponding to the other two matching locations are closer to
8 Lookup-Table Based Hyperspectral Data Compression Fig. 8.1 Previous image band. Co-located pixel has yellow background. Matching pixels have green background
Fig. 8.2 Current image band. Current pixel has yellow background. Pixels corresponding to the matching pixel have green backgrounds
173
336
335
314
335
314
335
319
327
316
315
317
315
328
315
325
319
322
334
329
314
329
324
317
315
328
339
323
339
328
332
331
335
335
324
325
327
320
332
327
335
330
350
339
324
333
325
333
325
Fig. 8.3 Lookup table
Index Value 314
328
315
332
316
335
317
333
the actual pixel value 325 than the nearest matching pixel value 332. This type of behavior of LUT method motivated the development of Locally Averaged Interband Scaling (LAIS)-LUT method [29], which uses a predictor to guide the selection between two LUTs. LAIS-LUT method works by first computing a LAIS estimate by scaling pixel co-located in the previous band. The LAIS scaling factor is an average of ratios between three neighboring causal pixels in the current and previous band: px;y 1;z px 1;y 1;z 1 px 1;y;z þ þ 3 px 1;y;z 1 px;y 1;z 1 px 1;y 1;z 1
(8.1)
LAIS scaling factor in (8.1) is used to compute an estimate for the current pixel: p
00 x;y;z
px;y 1;z px 1;y 1;z 1 px 1;y;z þ þ p ¼ 3 px 1;y;z 1 px;y 1;z 1 px 1;y 1;z 1 x;y;z
1
(8.2)
174 Fig. 8.4 LAIS estimates for LAIS-LUT
Fig. 8.5 Two lookup tables for LAIS-LUT
J. Mielikainen
Pixel LAIS Pixel Position Value Estimate (2,3)
324
320.1
(2,5)
327
321.9
(2,7)
332
316.2
index
1st LUT
2nd LUT
314
328
324
315
332
327
316
335
-
317
333
325
LAIS-LUT uses two LUTs, which are similar to the one used in the LUT method. The second LUT is updated with the past entries of the first LUT. The predictor returned by the LUT that is the closest one to the LAIS estimate is chosen as the predictor for the current pixel. If the LUTs return no match then the LAIS estimate is used as the estimated pixel value. We use the LUT example to illustrate the search process in LAIS-LUT. LAIS estimates for the three matching pixels in the previous example are shown in Fig. 8.4. Two LUTs corresponding to bands in Figs. 8.1 and 8.2 are shown in Fig. 8.5. Recall that the current pixel is p3,8,2 ¼ 325 and the causal pixels in the previous band are searched to find a match for the co-located pixel p3,8,1 ¼ 315. Out of the three matching pixels two are in LUTs (green background in Fig. 8.5). LAIS estimate (321.9) for 2nd LUT value 327 is closer than LAIS estimate (316.2) for the first LUT value 332. Therefore, pixel value from second LUT is used as the predictor for p’3,8,2 ¼ p2,5,2 ¼ 327.
8 Lookup-Table Based Hyperspectral Data Compression
175
5 Uniform Quantization of Co-Located Pixels In [24], a quantization of indices in LUT method was proposed. In LAIS-QLUT method a uniform quantization of the co-located pixels is performed before using them for indexing the LUTs. The use of quantization reduces the size of the LUTs by an order of magnitude A quantized interband predictor is formed by uniformly quantizing the colocated pixel px,y,z 1 before using it as an index to the LUT. Naturally, this reduces the size of the LUTs by the factor that is used in the uniform quantization. Except for a slightly simpler LAIS from [25] LAIS and an additional quantization step, LAIS-QLUT is the same algorithm as LAIS-LUT. The LAIS scaling factor in LAIS-QLUT is an average of ratios between three neighboring causal pixels in current and previous band: px 1;y;z þ px;y 1;z þ px 1;y 1;z 1 3 px 1;y;z 1 þ px;y 1;z 1 þ px 1;y 1;z 1
(8.3)
Thus, the corresponding LAIS estimate the current pixel is the following: p00 x;y;z ¼
px 1;y;z þ px;y 1;z þ px 1;y 1;z 1 p 3 px 1;y;z 1 þ px;y 1;z 1 þ px 1;y 1;z 1 x;y;z
(8.4)
1
LAIS in LAIS-QLUT requires a division operation and four addition operations compared to the three division, one multiplication, and two addition operations required by LAIS in LAIS-LUT. The search process in LAIS-QLUT will be illustrated using the same image bands are in the previous example. Quantized version of the previous image band is shown in Fig. 8.6 for a quantization factor 10. LAIS-Q estimates for two matching pixels are shown in Fig. 8.7. Two LUTs for LAIS-QLUT are shown in Fig. 8.8 for a quantization factor 10. The current pixel is p3,8,2 ¼ 325 and the causal pixels in the previous band are searched to find a match for quantized co-located pixel p3,8,1 / 10 ¼ 32. Two of matching pixels, which are in LUTs have LAIS-Q estimates of 328.2 for first LUT value 333 and 328.3 for second LUT value 325. The second LUT value is closer to the corresponding LAIS-Q estimate than the other one. Therefore, pixel value from the first LUT is used as the predictor for p’3,8,2 ¼ p3,6,2 ¼ 324.
Fig. 8.6 Quantized previous image band. Co-located pixel has yellow background. Matching pixels have green background
34
34
31
34
31
34
32
33
32
32
32
32
33
32
33
32
32
33
33
31
33
32
32
32
176 Fig. 8.7 LAIS estimates for LAIS-QLUT
Fig. 8.8 Two lookup table for LAIS-QLUT
J. Mielikainen
Pixel Position
Pixel Value
LAIS Estimate
(3,6)
325
328.3
(3,7)
333
328.2
index
1st LUT
2nd LUT
31
324
328
32
333
325
33
339
350
34
332
339
There are two separate variants of LAIS-QLUT. The first variant, The LAIS-QLUT-OPT method selects the optimal uniform quantization factor for each band. In order to find the optimal quantization factor, an exhaustive search of all possible quantization values is performed. Thus, the quantization factor selection is based on which quantization factor achieves the best compression efficiency for that specific band. The excessive time complexity of the LAIS-QLUT-OPT method could be decreased slightly by computing entropy of the residual image instead of actually encoding residuals for the determination of the optimal quantization factor. The second variant of LAIS-QLUT is called LAIS-QLUT-HEU and it uses constant quantization factors. The constant quantization factors are selected using a heuristic. The heuristic selects the constant quantization factors to be the bandwise mean values of the optimal quantization factors of an image set. A division operation required by the quantization represents the only increase in the time complexity of LAIS-QLUT-HEU compared to LAIS-LUT.
6 Multiband LUT In [26], LUT and LAIS-LUT method have been generalized to a multiband and multi-LUT method. In the extended method, the prediction of the current band relies on N previous bands. LUTs are defined on each of the previous bands
8 Lookup-Table Based Hyperspectral Data Compression
177
and each band contains M LUTs. Thus, there are NM different predictors to choose from. The decision among one of the possible prediction values is based on the closeness of the values contained in the LUTs to a reference prediction. Two different types of purely spectral multiband prediction estimates were proposed for. One of the reference predictors is crisp and the other one is fuzzy. The first method is S-RLP [15]. The method partitions image bands into blocks, and a predictor, out of a set of predictors, is selected for prediction. In the S-FMP method [15] the causal neighborhoods of each pixel are clustered using fuzzy-cmeans clustering. For each of the clusters, an optimal linear predictor is computed from the values, the membership degrees of which exceed a threshold. The final estimate is computed as a weighted sum of the predictors, where the weights are the membership degrees. The LUT based compression methods based on S-RLP and S-FMP are denoted as S-RLP-LUT and S-FMP-LUT, respectively.
7 Experimental Results Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) is an airborne hyperspectral system collecting spectral radiance in 224 contiguous spectral bands with wavelengths from 370 to 2,500 nm. The AVIRIS instrument consists of four spectrometers that view a 20-m2 spot on the ground from a flight altitude of 20 km. This spot is simultaneously viewed in all the spectral bands. A spatial image is formed by moving the spectrometers perpendicular to the direction of the aircraft [27]. Experimental results are shown for two different AVIRIS data sets. The first data set consists of four calibrated radiance images from 1997 AVIRIS sample data product. The AVIRIS images are from the following four different areas: Cuprite, NV; Jasper Ridge, CA; Lunar Lake, NV; and Moffett Field, CA. They are the most widely used data for benchmarking hyperspectral image compression algorithms. Image features and the number of lines are listed in Table 8.1. Each image contains 614 samples/line and they are stored as 16-bit signed integers. A gray scale image of Moffett Field image can be seen in Fig. 8.9. Newer data set was acquired on 2006. A new AVIRIS data set consists of five calibrated and uncalibrated 16-bit images from Yellowstone, WY and two 12-bit uncalibrated images one from Hawaii and one from Maine. Summary of the new Consultative Committee for Space Data Systems (CCSDS) AVIRIS data is given in Table 8.2. Each image is a 512-line scene containing 224 spectral bands. An example of a scene can be seen in Fig. 8.10 in the form of a false color image of calibrated Yellowstone scene 11. Table 8.1 The standard 1997 AVIRIS images [11] Site Features Cuprite Geological features Jasper Ridge Vegetation Lunar Lake Calibration Moffett Field Vegetation, urbar, water
Lines 2,206 2,586 1,431 2,031
178
J. Mielikainen
Fig. 8.9 Gray scale image of Moffett Field image from AVIRIS 1997 image set
This AVIRIS data is a part of the CCSDS data set, which is used to assess the performance of hyperspectral compression algorithms. Table 8.3 shows results for the NN method. The first column depicts the length of the search window; 0 lines means that only the current line is searched. The following columns are bit rates in bits/pixel for the four test images and the average, respectively. When the search window’s length is equal to the length of image, the method naturally predicts the same values as the LUT method. These results show
8 Lookup-Table Based Hyperspectral Data Compression Table 8.2 AVIRIS images included in the CCSDS test set [11] Site Scene numbers Year Samples/line Yellowstone 0,3,10,11,18 2006 677 Yellowstone 0,3,10,11,18 2006 680 Hawaii 1 2001 614 Maine 10 2003 680
179
Bit depth 16 16 12 12
Type Calibrated Uncalibrated Uncalibrated Uncalibrated
Fig. 8.10 False color image of calibrated Yellow stone 11 from CCSDS AVIRIS data set
Table 8.3 Compression results in bits/pixel for calibrated AVIRIS 1997 test images in bits per pixel # of lines Cuprite Jasper ridge Lunar lake Moffett field 0 5.69 5.84 5.78 6.02 1 5.41 5.63 5.50 5.80 2 5.29 5.50 5.33 5.65 4 5.05 5.35 5.14 5.48 8 4.89 5.21 4.98 5.32 16 4.79 5.10 4.88 5.21 32 4.72 5.03 4.79 5.14 64 4.69 5.00 4.75 5.10 128 4.68 4.98 4.73 5.08 256 4.66 4.97 4.72 5.06 512 4.66 4.97 4.72 5.05 1,024 4.65 4.95 4.71 5.05
180
J. Mielikainen
Table 8.4 Compression results in bits/pixel for calibrated AVIRIS 1997 test images in bits per pixel Cuprite Jasper ridge Lunar lake Moffett field Average JPEG-LS 7.66 8.38 7.48 8.04 7.89 Diff. JPEG-LS 5.50 5.69 5.46 5.63 5.57 3D-CALIC 5.23/5.39 5.19/5.37 5.18/5.32 4.92/5.05 5.19/5.28 BH –/5.11 –/5.23 –/5.11 –/5.26 –/5.18 M-CALIC 4.97/5.10 5.05/5.23 4.88/5.02 4.72/4.89 4.98/5.06 SLSQ-OPT 4.94/5.08 4.95/5.08 4.95/5.08 4.98/5.10 4.96/5.09 CCAP –/4.92 –/4.95 –/4.97 – – KSP –/4.88 –/4.95 –/4.89 –/4.92 –/4.91 FL# 4.82 4.87 4.82 4.93 4.86 NPHI 4.79 4.89 4.97 4.79 4.86 C-DPCM –/4.68 –/4.62 –/4.75 –/4.62 –/4.67 S-RLP 4.69 4.65 4.69 4.67 4.67 S-FMP 4.66 4.63 4.66 4.63 4.64 LUT 4.66 4.95 4.71 5.05 4.84 LAIS-LUT 4.47 4.68 4.53 4.76 4.61 LAIS-QLUT-HEU 4.30 4.62 4.36 4.64 4.48 LAIS-QLUT-OPT 4.29 4.61 4.34 4.63 4.47 S-RLP-LUT 3.92 4.05 3.95 4.09 4.00 S-FMP-LUT 3.89 4.03 3.92 4.05 3.97 IP3-BPS 3.76 4.06 3.79 4.06 3.92
that limiting the search window size significantly affects the performance of the NN method compared to the full search. Thus, a large search window is necessary in order to achieve good compression ratios. Table 8.4 shows compression results for AVIRIS 1997 data. The results are reported for the band-interleaved-by-line (BIL) and band-sequential (BSQ) formats. In the BIL format, the current line, along with the two previous lines, is available. For BSQ data, the current band and several previous bands are available for processing. The LUT family does not benefit from the BSQ data format. This is due to two factors. First, LUT and LAIS-LUT methods only utilize one previous band. Second, LAIS-LUT methods need only the data from the current and previous image lines. Those lines were already provided by the BIL data format. Most compression method exhibit identical compression results for both BIL and BSQ data. Only one bits/pixel value is shown for those methods. For the other methods both BIL and BSQ results are provided. The results for the two different data formats are separated by a forward-slash and dash denotes unavailable results. Differential JPEG-LS computes the difference between each band and the previous band before running JPEG-LS on residual data. Experimental results show that LUT based algorithms work extremely well for calibrated AVIRIS 1997 data. Even the low time complexity LAIS-LUT and QLAIS-LUT variants have close to the state-of-the-art compression ratios. IP3BPS method takes ten times longer than LUT and five times longer than LAIS-LUT or LAIS-QLUT-HEU to compress AVIRIS image [13].
8 Lookup-Table Based Hyperspectral Data Compression
181
Table 8.5 Compression results in bits/pixel for 16-bit raw CCSDS AVIRIS test images in bits per pixel Algorithm Scene 0 Scene 3 Scene 10 Scene 11 Scene 18 Average JPEG-LS 9.18 8.87 7.32 8.50 9.30 8.63 BG 6.46 6.31 5.65 6.05 6.40 6.17 A1 6.92 6.78 6.10 6.53 6.92 6.65 LUT 7.13 6.91 6.25 6.69 7.20 6.84 LAIS-LUT 6.78 6.60 6.00 6.30 6.82 6.50 FL# 6.20 6.07 5.60 5.81 6.26 5.99 IP3 6.20 6.08 5.56 5.81 6.25 5.98 C-DPCM-20 5.88 5.71 5.20 5.52 5.75 5.61 C-DPCM-80 5.82 5.65 5.17 5.47 5.69 5.56
Table 8.6 Compression results in bits/pixel for 12-bit raw CCSDS AVIRIS test images in bits per pixel Algorithm Hawaii Maine Average JPEG-LS 4.58 4.50 4.54 A1 3.49 3.65 3.57 LUT 3.27 3.44 3.36 LAIS-LUT 3.05 3.19 3.12 BG 3.03 3.17 3.10 IP3 2.55 2.68 2.62 FL# 2.58 2.63 2.61 C-DPCM-20 2.43 2.57 2.50 C-DPCM-80 2.38 2.52 2.45
The LUT method requires a full LUT for each band. Assuming 16-bit LUTs, each LUT’s memory requirements are roughly equivalent to 107 lines of an AVIRIS image data. The LUT’s memory requirements are independent of the spatial size of the image. Therefore, the relative size of the LUTs compared to the image gets smaller as the spatial size of the image gets larger. For our test images, the amount of the memory required by LUTs is 4–7% of the memory used by the image. The average quantization factor for LAIS-QLUT-HEU was 28. Thus, the average LUT memory requirement is roughly equivalent to four lines of AVIRIS image data compared to 107 lines of data in the original LUT method. We have also experimented with the optimization of the quantization factors for each image instead of for each band. That procedure gave a quantization factor of ten for all the test images. The average bit rate was 4.60 bits/pixel. This compares unfavorably to the 4.47 bits/pixel average bit rate of LAIS-QLUT-HEU. Therefore, separate bandwise quantization factors are worthwhile. Tables 8.5–8.7 depict compression results for new AVIRIS data in bits per pixel for various different compression methods. C-DPCM-20 and C-DPCM-80 refer to the prediction length 20 and 80 for C-DPCM, respectively. A modified C-DPCM method uniformly quantizes coefficients to 12 bits instead of 16 bits in the original C-DPCM.
182
J. Mielikainen
Table 8.7 Compression results in bits/pixel for calibrated CCSDS AVIRIS test images in bits per pixel Algorithm Scene 0 Scene 3 Scene 10 Scene 11 Scene 18 Average JPEG-LS 6.95 6.68 5.19 6.24 7.02 6.42 A1 4.81 4.69 4.01 4.41 4.77 4.54 LUT 4.81 4.62 3.95 4.34 4.84 4.51 LAIS-LUT 4.48 4.31 3.71 4.02 4.48 4.20 BG 4.29 4.16 3.49 3.90 4.23 4.01 FL# 3.91 3.79 3.37 3.59 3.90 3.71 IP3 3.81 3.66 3.13 3.45 3.75 3.56 C-DPCM-20 3.61 3.43 2.97 3.28 3.49 3.36 C-DPCM-80 3.53 3.36 2.93 3.22 3.43 3.29
The results for uncalibrated CCSDS AVIRIS test data in Tables 8.5 and 8.6 show that LUT-based methods lose their performance advantage when applied to uncalibrated data. Moreover, the results in Table 8.7 show that LUT-based algorithms that exploit calibration artifacts in AVIRIS 1997 images have no performance advantage on the calibrated CCSDS AVIRIS images.
8 Conclusions An overview of the lookup table (LUT) based lossless compression methods for hyperspectral images have been presented in this chapter. Experimental results on AVIRIS data showed that the LUT based algorithms work extremely well for old calibrated AVIRIS data. Even the low-complexity LAIS-LUT and QLAIS-LUT variants have close to the state-of-the-art compression ratios. LUT-based methods exploit artificial regularities that are introduced by the conversion of raw data values to radiance units [11]. The calibration-induced artifacts are not present in the newer AVIRIS images in Consultative Committee for Space Data Systems (CCSDS) test set. Thus, LUT based method do not work as well on raw or the newer AVIRIS images in 2006, which use new calibration measures. Acknowledgement This work was supported by the Academy of Finland.
References 1. A. Bilgin, G. Zweig, and M. Marcellin, “Three-dimensional image compression with integer wavelet transforms,” Appl. Opt., vol. 39, no. 11, pp. 1799–1814, Apr. 2000. 2. B. Baizert, M. Pickering, and M. Ryan, “Compression of hyperspectral data by spatial/spectral discrete cosine transform,” in Proc. Int. Geosci. Remote Sens. Symp., 2001, vol. 4, pp. 1859–1861, doi: 10.1109/IGARSS.2001.977096.
8 Lookup-Table Based Hyperspectral Data Compression
183
3. J. Mielikainen and A. Kaarna, “Improved back end for integer PCA and wavelet transforms for lossless compression of multispectral images,” in Proc. 16th Int. Conf. Pattern Recog., Quebec City, QC, Canada, 2002, pp. 257–260, doi: 10.1109/ICPR.2002.1048287. 4. M. Ryan and J. Arnold, “The lossless compression of AVIRIS images vector quantization,” IEEE Trans. Geosci. Remote Sens., vol. 35, no. 3, pp. 546–550, May 1997, doi: 10.1109/ 36.581964. 5. J. Mielikainen and P. Toivanen, “Improved vector quantization for lossless compression of AVIRIS images,” in Proc. XI Eur. Signal Process. Conf., Toulouse, France, Sep. 2002, pp. 495–497. 6. G. Motta, F. Rizzo, and J. Storer, “Partitioned vector quantization application to lossless compression of hyperspectral images,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., Jul. 2003, vol. 1, pp. 553–556, doi: 10.1109/ICME.2003.1220977. 7. S. Tate, “Band ordering in lossless compression of multispectral images,” IEEE Trans. Comput., vol. 46, no. 4, pp. 477–483, Apr. 1997, doi: 10.1109/12.588062. 8. P. Toivanen, O. Kubasova, and J. Mielikainen, “Correlation-based bandordering heuristic for lossless compression of hyperspectral sounder data,” IEEE Geosci. Remote Sens. Lett., vol. 2, no. 1, pp. 50–54, Jan. 2005, doi: 10.1109/LGRS.2004.838410. 9. J. Zhang and G. Liu, “An efficient reordering prediction based lossless compression algorithm for hyperspectral images,” IEEE Geosci. Remote Sens. Lett., vol. 4, no. 2, pp. 283–287, Apr. 2007, doi: 10.1109/LGRS.2007.890546. 10. A. Abrardo, M. Barni, E. Magli, F. Nencini, “Error-Resilient and Low-Complexity On-board Lossless Compression of Hyperspectral Images by Means of Distributed Source Coding,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 4, pp. 1892–1904, 2010, doi:10.1109/ TGRS.2009.2033470. 11. A. B. Kiely, M. A. Klimesh, “Exploiting Calibration-Induced Artifacts in Lossless Compression of Hyperspectral Imagery,” IEEE Trans. Geosci. Remote Sens., vol. 47, no. 8, pp. 2672–2678, 2009, doi:10.1109/TGRS.2009.2015291. 12. M. Slyz, L. Zhang, “A block-based inter-band lossless hyperspectral image compressor,” in Proc. of IEEE Data Compression Conference, pp. 427–436, 2005, doi: 10.1109/DCC.2005.1. 13. C.-C. Lin, Y.-T. Hwang., “An Efficient Lossless Compression Scheme for Hyperspectral Images Using Two-Stage Prediction”, vol. 7, no. 3, pp. 558–562, 2010, doi:10.1109/ LGRS.2010.2041630. 14. J. Mielikainen, P. Toivanen, “Clustered DPCM for the Lossless Compression of Hyperspectral Images”, IEEE Trans. Geosci. Remote Sens., vol. 41, no. 12, pp. 2943–2946, 2003 doi:10.1109/TGRS.2003.820885. 15. B. Aiazzi, L. Alparone, S. Baronti, and C. Lastri, “Crisp and fuzzy adaptive spectral predictions for lossless and near-lossless compression of hyperspectral imagery,” IEEE Geosci. Remote Sens. Lett., vol. 4, no. 4, pp. 532–536, Oct. 2007, 10.1109/LGRS.2007.900695. 16. E. Magli, G. Olmo, and E. Quacchio, “Optimized onboard lossless and near-lossless compression of hyperspectral data using CALIC,” IEEE Geosci. Remote Sens. Lett., vol. 1, no. 1, pp. 21–25, Jan. 2004, doi:10.1109/LGRS.2003.822312. 17. F. Rizzo, B. Carpentieri, G. Motta, and J. Storer, “Low-complexity lossless compression o hyperspectral imagery via linear prediction,” IEEE Signal Process. Lett., vol. 12, no. 2, pp. 138–141, Feb. 2005, doi:10.1109/LSP.2004.840907. 18. J. Mielikainen and P. Toivanen, “Parallel implementation of linear prediction model for lossless compression of hyperspectral airborne visible infrared imaging spectrometer images,” J. Electron. Imaging, vol. 14, no. 1, pp. 013010-1–013010-7, Jan.–Mar. 2005, doi:10.1117/ 1.1867998. 19. H. Wang, S. Babacan, and K. Sayood, “Lossless Hyperspectral-Image Compression Using Context-Based Conditional Average,” IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 12, pp. 4187–8193, Dec. 2007, doi:0.1109/TGRS.2007.906085.
184
J. Mielikainen
20. M. Slyz and L. Zhang, “A block-based inter-band lossless hyperspectral image compressor,” in Proc. Data Compression Conf., Snowbird, UT, 2005, pp. 427–436, doi:10.1109/ DCC.2005.1. 21. S. Jain and D. Adjeroh, “Edge-based prediction for lossless compression of hyperspectral images,” in Proc. Data Compression Conf., Snowbird, UT, 2007, pp. 153–162, doi:10.1109/ DCC.2007.36. 22. J. Mielikainen, “Lossless compression of hyperspectral images using lookup tables,” IEEE Sig. Proc. Lett., vol. 13, no. 3, pp. 157–160, 2006, doi:10.1109/LSP.2005.862604. 23. E. Magli, “Multiband lossless compression of hyperspectral images,” IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 4, pp. 1168–1178, Apr. 2009, doi:10.1109/ TGRS.2008.2009316. 24. J. Mielikainen, P. Toivanen, “Lossless Compression of Hyperspectral Images Using a Quantized Index to Lookup Tables,” vol. 5, no. 3, pp. 474–477, doi:10.1109/LGRS.2008.917598. 25. J. Mielikainen, P. Toivanen, and A. Kaarna, “Linear prediction in lossless compression of hyperspectral images,” Opt. Eng., vol. 42, no. 4, pp. 1013–1017, Apr. 2003, doi:10.1117/ 1.1557174. 26. B. Aiazzi, S. Baronti, S., L. Alparone, “Lossless Compression of Hyperspectral Images Using Multiband Lookup Tables,” IEEE Signal Processing Letters, vol. 16, no. 6, pp. 481–484. Jun. 2009, doi:10.1109/LSP.2009.2016834, 0.1109/LSP.2009.2016834. 27. W. Porter and H. Enmark, “A system overview of the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS),” Proc. SPIE, vol. 834, pp. 22–31, 1997 28. X. Wu and N. Memon, “Context-based lossless interband compression—Extending CALIC,” IEEE Trans. Image Process., vol. 9, no. 6, pp. 994–1001, Jun. 2000, doi:10.1109/83.846242. 29. B. Huang, Y. Sriraja, “Lossless compression of hyperspectral imagery via lookup tables with predictor selection,” in Proc. SPIE, vol. 6365, pp. 63650L.1–63650L.8, 2006, doi:10.1117/ 12.690659.
Chapter 9
Multiplierless Reversible Integer TDLT/KLT for Lossy-to-Lossless Hyperspectral Image Compression Jiaji Wu, Lei Wang, Yong Fang, and L.C. Jiao
1 Introduction Hyperspectral images have wide applications nowadays such as in atmospheric detection, remote sensing and military affairs. However, the volume of a hyperspectral image is so large that a 16bit AVIRIS image with a size 512 512 224 will occupy 112 M bytes. Therefore, efficient compression algorithms are required to reduce the cost of storage or bandwidth. Lossy-to-lossless compression will be of great importance in telemedicine and satellite communications for legal reasons and research requirements. To realize scalable coding, most of the state-of-the-art compression methods adopt three dimensional discrete wavelet transform (3D-DWT) [1–3] or wavelet transform/ karhunen-loeve transform (DWT/KLT) [4–6], where a 9/7 floating-point filter (9/7F filter) is always used for de-correlation in lossy compression. Lossless compression schemes include methods based on vector quantization (VQ), prediction, integer transforms and so on. Although prediction-based methods perform well, they do not have the ability to perform progressive lossy-to-lossless compression since this depends on transform-based methods [7]. Sweldens [8, 9] proposed a lifting scheme for the realization of wavelet transforms. Bilgin et al. [10] introduced a reversible integer wavelet transform method for 3D image compression. Xiong et al. [11] applied 3D integer wavelet transforms for medical image compression and pointed out that the transform has to be unitary to achieve good lossy coding performance. Some researchers have studied integer KLT for spectral decorrelation. Hao et al. [12] proposed reversible integer KLT (RKLT) and Galli et al. [13] improved it. However, in the spatial domain, integer wavelet transforms are still commonly
J. Wu (*) • L. Wang • L.C. Jiao School of Electronic Engineering, Xidian University, xi’an, China e-mail: [emailprotected] Y. Fang College of Information Engineering, Northwest A&F University, yangling, China B. Huang (ed.), Satellite Data Compression, DOI 10.1007/978-1-4614-1183-3_9, # Springer Science+Business Media, LLC 2011
185
186
J. Wu et al.
B-F 0
1 B-F
2
B-F 3
4
5
B-F
Fig. 9.1 B-F standards for basis functions of LOT; one B-F imposes on three neighboring segments but just output one segment of values
used. A drawback of the wavelet-based compression method is that 53DWT is usually applied instead of 97DWT in lossy-to-lossless compression schemes, and this will lead to performance degradation. Another disadvantage of DWT is that it cannot compete with DCT, due to the constraint of CPU performance and computer memory, especially in real-time and low-complexity applications, because the computing complexity of DWT increases exponentially with image size [14]. The computational complexity of DWT increases when the image size increases because of the global transform [14]. However, DCT has its own special advantages such as low memory cost, flexibility at block by block level and parallel processing. DCT is approximately equal to the KLT basis matrix for the first-order Markov process while image segments always satisfy this condition. Therefore, DCT performs well in image decorrelation and it is widely adopted by image/video compression standards such as JPEG, MPEG, and H.26X. Although the DCTbased coding method has been a popular method for image and video compression, a key problem of this type of coding at low bit rates is the so-called “block effect”. The reason is because DCT-based coding always independently processes each block. In order to reduce the block effect of DCT compression, some deblocking methods based on filtering were proposed [15, 16], in which some low-pass filters were applied to the boundary pixels. However, filtering-based methods usually blur image content. To resolve the problem of the block effect of DCT, lapped orthogonal transform (LOT) was proposed by Cassereau [17] and an analytical solution has been given by Malvar [18]. LOT improves DCT by designing basis functions which impose themselves on neighboring segments, as depicted in Fig. 9.1. Correlations between neighboring segments can be explored in this way, and discontinuity between reconstructed segments can be reduced [19]. In lossy compression, although LOT can reduce block effects efficiently, the lapped filtering of LOT has to follow behind the DCT. For this reason, forward LOT is difficult to make compatible with a DCTbased coding standard. To overcome the disadvantage of traditional LOT, Tran et al. [20, 21] have designed a family of time domain lapped transforms (TDLT) by adding various preand post-filters in the existing block-based architecture in the time domain. Tran’s algorithm can achieve competitive compression performance compared with DWTbased methods while reducing or even eliminating the block artifacts to guarantee good visual quality. TDLT, a combination of pre- and post-filters with DCT transform, can be illustrated in this way; the inputs of DCT and the outputs of the
9 Multiplierless Reversible Integer TDLT/KLT. . .
187
inverse DCT are processed by a pre-filter and post-filter, respectively. The pre-filter for TLDT is placed in front of the forward DCT, so the TDLT is easily made compatible with current DCT-based coding standards. The filtering process is conducted on the two neighboring blocking coefficients. The function of the prefilter is to make the input data of each DCT block as homogenous as possible, like a flattening operation, whereas the function of the post-filter is to reduce the block artifacts. In [20], lossy compression has been realized by using different versions of TDLT based on various decompositions of the filtering matrix, and lossless compression has also been realized by using a reversible transform for lifting-based filters and multiplier-less approximations of DCT, known as binDCT in [22]. BinDCT is realized by quantizing the transform coefficients of conventional plane rotation-based factorizations of the DCT matrix, and can be implemented using only binary shift and addition operations. Microsoft has developed a TDLT-based image coding technique called HDPhoto [23, 24], which enables reversible compression by using lifting scheme. HD-Photo has been taken as the basis technique for JPEG-XR, which is a new compression format supporting high dynamic range and promising significant improvements in image quality and performance for end-to-end digital photography. In HD-Photo, a hierarchical lapped biorthogonal transform (LBT) is adopted and the Huffman coding is performed in chunks organized as a function of resolution [23]. Both lossy and lossless compression can be realized by HDPhoto and this is one of its advantages over JPEG, which needs two coders in different applications. In HD-Photo, the LBT is realized by factorizing the core transform and overlap filtering into rotation operators and is implemented using a lifting structure to promise a reversible integer transform [25]. In addition, the new compression scheme retains advantages such as in-place calculation, amenability to parallelized implementation, flexibility and adaptivity on the block level and so on. Although TDLT performs even better than DWT does in energy compatibility and lossy compression, it does not perform well in the lossless compression where the reversible transform is required. In fact, for hyperspectral image compression, a completely reversible transform method is often required to realize lossy-to-lossless coding. In this chapter, we take a practical and innovative approach to replace integer DWT with integer reversible time domain lapped transform (RTDLT) in the spatial domain, and RKLT is applied in the spectral dimension. Here this RTDLT and RKLT are realized by an improved matrix factorization method. RTDLT can realize integer reversible transform and hence we have adopted a progressive lossy-to-lossless hyperspectral image compression method based on RTDLT and RKLT. Block transforming coefficients in the spatial domain are reorganized into sub-band structures so as to be coded by wavelet-based coding methods. In addition, an improved 3D embedded zero-block coding method used to code transformed coefficients is integrated in this work. Moverover, we also extend RTDLT to 3D reversible integer lapped transform (3DRLT), which can replace 3D integer WT and realize progressive lossy- to-lossless
188
J. Wu et al.
compression and performs better than 3D-WT in most cases. 3D-RLT is implemented at block-level and has a fixed transform basis matrix. It is therefore suitable for memory-limited systems or on board spacecraft where component constraints are significant. Our proposed methods retain the character of scalability in reconstruction quality and spatial resolution so, at the decoder, observers can review the whole image from inferior quality to the completely reconstructed image. Experimental results show that the proposed methods perform well in both lossy and lossless compression. To reduce complexity, our proposed methods are implemented using shift and addition without any multiplier, with the help of multi-lifting.
2 Multi-lifting Scheme We adopt a lifting scheme to realize the reversible integer transform since lifting schemes have been widely used and have many advantages such as (a) fast implementation, (b) in-place calculations, (c) immediate access to inverse transform, (d) a natural understanding of original complex transforms. Firstly, we give a simple review of our lifting scheme.
2.1
Lifting Scheme and Applications
To realize a reversible integer transform, traditional algorithms always adopt duallifting, as depicted in Fig. 9.2 where xi and yi stand for input and output signal, respectively and p and ustand for prediction and update coefficients, respectively. For instance, Daubechies and Sweldens proposed a lifting-based wavelet transform method [9], Chen et al. [26] proposed Integer DCT using the Walsh-Hadamard Transform (WHT) and lifting scheme, Abhayaratne [27] proposed an N-point integer-to-integer DCT (I2I-DCT) by applying recursive methods and lifting techniques, where N is power of 2 and Liang et al. [22] proposed two kinds of fast multiplier-less approximations of the DCT called binDCT, also with the lifting scheme. In the HD-Photo process from Microsoft, the LBT is realized by factorizing the core transform and overlap filtering into rotational operators, as depicted in Fig. 9.3; the rotation operations are also realized by a dual-lifting structure [25]. Matrix lifting has also been proposed by Li [28], andis applied in embedded audio codec. Cheng et al. have introduced the properties of the lifting matrix, based on which a new family of lapped biorthogonal transforms has been designed [29]. Lifting-based transforms can realize completely reversible integer transforms and hence can be applied in lossless compression. At the same time, lifting is lower in CPU cost and memory usage since it allows a fully in-place calculation.
9 Multiplierless Reversible Integer TDLT/KLT. . .
189
Fig. 9.2 Dual lifting structure: (a) Forward. (b) Inverse
Fig. 9.3 Transform structure of LBT in HD-Photo
2.2
Reversible Integer Transform based on Multi-lifting Scheme
In this chapter, we adopt a multi-lifting scheme, which is an extension of traditional dual lifting. Multi-lifting is the same concept as matrix lifting, which has been used in the JPEG2000 color transform. If a 2 2 matrix U is an upper triangular matrix with the diagonal elements equal to 1, then Y ¼ UX can be realized as (9.1) from integer to integer; it can also be implemented by the lifting depicted in Fig. 9.4a. In addition, X ¼ U 1Y can be realized using inverse-lifting as depicted in Fig.9.4b.
y0 y1
¼
1 0
p 1
x0 x1
!
y0 ¼ x0 þ bpx1 c x1 ¼ y1 ; y1 ¼ x1 x0 ¼ y0
bpy1 c
(9.1)
In our proposed transforming scheme, 4-order, 8-order and 16-order lifting are used. For example, the 4-point reversible integer-to-integer transforms Y ¼ UX and its inverse transform can be realized as (9.2), where U is an upper triangular
190
J. Wu et al.
Fig. 9.4 Reversible integer to integer transform based on lifting: (a) Forward integer lifting. (b) Inverse integer lifting
Fig. 9.5 Forward and inverse multi-lifting schemes
matrix of size 4. The multi-lifting implementation is also depicted in Fig. 9.5. What should be noted isP that the same rounding operation is implied on the sum of a set of N 1 multiplications ( j¼iþ1 ui;j xj ) in the forward and inverse transforms to guarantee completely reversible integer-to-integer transforms. y0 ¼ x0 þ bu01 x1 þ u02 x2 þ u03 x3 c x3 ¼ y3 ; x2 ¼ y2 bu23 x3 c; y1 ¼ x1 þ bu12 x2 þ u13 x3 c x1 ¼ y1 bu12 x2 þ u13 x3 c; y2 ¼ x2 þ bu23 x3 c x0 ¼ y0 bu01 x1 þ u02 x2 þ u03 x3 c; y3 ¼ x3
(9.2)
9 Multiplierless Reversible Integer TDLT/KLT. . .
191
3 Reversible Integer to Integer TDLT/RKLT In this section, we will introduce the implementation of RTDLT/RKLT which is based on the matrix factorization method [30]. We will first introduce techniques of integer-to-integer transform by multi-lifting based on matrix factorization. We will then demonstrate details in the design of RTDLT/RKLT.
3.1
Matrix Factorization and Multi-lifting
Based on matrix factorization theory [30], a nonsingular matrix can be factorized into a product of at most three triangular elementary reversible matrices (TERMs). If the diagonal elements of the TERM are equal to 1, the reversible integer to integer transform can be realized by multi-lifting. First, we review how to realize an approximation of floating-point to integer for a transform basis using matrix factorization theory. Suppose A is a transform basis matrix, a1;1
ð0Þ
a1;2
B B ð0Þ B a2;1 B A¼B B .. B . @
ð0Þ a2;2
ð0Þ aN;1
ð0Þ
a1;N
...
ð0Þ a2;N
..
.. . ð0Þ aN;2
ð0Þ
...
1
C C C C C: .. C . C A
.
(9.3)
ð0Þ
. . . aN;N
There exists a permutation matrix P1 which can make P11;N not equal to zero: 0
ð1Þ
p1;1
B B ð1Þ B p2;1 B P1 A ¼ B B .. B . @ ð1Þ
pN;1
p1;2
ð1Þ
p2;2
ð1Þ
.. .
..
.
ð1Þ
ð1Þ
p1;N
1
C ð1Þ C p2;N C C C: .. C . C A
(9.4)
ð1Þ
pN;2
pN;N
There is an operator s1 satisfying the following formula, ð1Þ
P1;1
ð1Þ
s1 P1;N ¼ 1:
(9.5)
Formula (9.5) can be rewritten as, ð1Þ s1 ¼ P1;1
ð1Þ 1 =P1;N :
(9.6)
192
J. Wu et al.
Then
2
P1 AS1 ¼ P1 A4
1 s1
I 0
2
1
6 6 ð1Þ 6 p2;1 5¼6 6 6 .. 1 6. 4 3
ð1Þ
s1 p2;N
ð1Þ
ð1Þ
pN;1
2
ð1Þ
p2;2
ð1Þ
.. .
.. ð1Þ
s1 pN;N
1 6 ð1Þ ð1Þ 6 s1 p2;N p2;1 6 L1 ¼ 6 6 .. 6 . 4 ð1Þ ð1Þ s1 pN;N pN;1
p1;2
1 I
1
3
7 ð1Þ 7 p2;N 7 7 7: 7 .. 7 . 5 ð1Þ pN;N
.
3
pN;2
.. .
ð1Þ
p1;N
7 7 7 7: 7 7 5
(9.7)
(9.8)
A Gaussian matrix which satisfies formula (9.8) will produce the following result, 3 3 2 2 ð2Þ ð2Þ 1 1 a1;2 a1;N 7 7 6 6 s pð1Þ pð1Þ 1 7 6 6 1 2;N ð2Þ ð2Þ 7 2;1 7 6 0 a2;2 6 a2;N 7 L1 P1 AS1 ¼ 6 7P1 AS1 ¼ 6 7: .. .. 7 6 6 7 . I . 5 5 4 4 ð2Þ ð2Þ ð1Þ ð1Þ 0 aN;2 aN;N s1 pN;N pN;1 1 (9.9)
Denoting Að1Þ ¼ L1 P1 AS1 , according to the method described above, there exists ð2Þ a P2 which satisfies P2;N 6¼ 0, together with an s2 which satisfies formula (9.10). ð2Þ
ð2Þ
s2 P2;N ¼ 1:
P2;2
(9.10)
Following this recursive process, Pk ; sk ; Lk ðk ¼ 1; 2; NÞ can be determined. Accordingly, we arrive at the following formula:
LN 1 PN
1
L2 P2 L1 P1 AS1 S2 SN
2
1
6 6 6 ¼6 6 4
1 0
a1;2
ðN 1Þ
a1;N
1
a2;N
ðN 1Þ ðN 1Þ
ðN 1Þ
aN;N
3
7 7 7 7 ¼ DR U; 7 5
(9.11)
where DR ¼ diagð1; 1; ; 1; eiy Þ:
(9.12)
9 Multiplierless Reversible Integer TDLT/KLT. . .
2
ðN 1Þ
a1;2
1
193
ðN 1Þ
a1;N
3
7 6 6 ðN 1Þ 7 1 a2;N 7 6 0 U¼6 7: 6 7 5 4 0 0 1
(9.13)
Setting L
1
¼ LN 1 ðPN 1 LN 2 PTN 1 Þ ðPN 1 PN P T ¼ PN 1 PN S
1
2
2
P2 L1 PT2 PT3 PTN 2 PTN 1 Þ: (9.14)
P2 P1 :
(9.15)
¼ S1 S2 SN 1 :
(9.16)
We conclude that L 1 PT AS
1
¼ DR U:
(9.17)
We then get A ¼ PLUS. What should be noted is that the factorization is not unique and different factorizations affect the error between the integer approximation transform and the original transform. This will also affect the intrinsic energy-compacting capability of the original transform, so that the error should be reduced as much as possible. Quasi-complete pivoting is suggested in the progress of matrix factorization [13]. In our experiments, we found this method to be very effective in reducing the error and enhancing the stability, so achieving better integer approximation to the original floating-point transform. The improved implementation proceeds as follows: Create a new matrix Sc1 using the transform basis matrix A: 2
6 6 Sc1 ¼ 6 4
ð0Þ
ð0Þ
ða1;1 Þ
1=a1;2
..
.. . ð0Þ
ðaN;1 Þ
ð0Þ
1=aN;2
ð0Þ
1=a1;N .. .
.
ð0Þ
ða1;1 Þ ð0Þ
ð0Þ
ðaN;1 Þ
1=aN;N
3
7 7 7: 5
(9.18)
Choose a parameter s1 which has a minimum absolute value in matrix Sc1 . This method is different from the traditional method which places a nonzero element at the end of every row of the transform basis matrix A in order to perform the calculation. ð0Þ
s1 ¼ minfSc1 g ¼ ðai;1
ð0Þ
1Þ=ai;j :
(9.19)
194
J. Wu et al.
If i is not equal to 1, that is, s1 is not located in the first row of A, then the i-th row must be permuted to the first row using a permutation matrix P1 . 0
ð1Þ
ð1Þ
q1;1
q1;2
B B qð1Þ B 2;1 P1 A ¼ B B . B . @ .
q2;2
ð1Þ
.. .
..
ð1Þ qN;1
ð1Þ qN;2
.
ð1Þ
q1;N
1
C ð1Þ q2;N C C C: .. C C . A
(9.20)
ð1Þ
qN;N
Create a matrix S1 which has the following shape: 2
1 0 0 0
60 6 6 6. 6 .. 6 S1 ¼ 6 6 s1 6 6 . 6 . 4 .
3
1 0 07 7 7 .. 7 .. . . .7 . . 7 7: 0 1 07 7 . . .. 7 .. 7 . .5 .
(9.21)
0 0 0 1
Then according to formula (9.18) to (9.21), we can get a more stable A ¼ PLUS. It can be proved that, based on the theory described above, matrix A has a unit TERM factorization of A ¼ PLUS if and only if det A ¼ det P ¼ 1, where L and S are unit lower triangular matrices, U is an unit upper triangular matrix and P is a permutation matrix. The type-II DCT [31] in one dimension is given by the following formula: XC ðkÞ ¼ ek
rffiffiffiffi N 1 2X N
n¼0
xðnÞ cosðð2n þ 1Þ
pk Þ; 2N
(9.22)
pffiffiffi for k ¼ 0. . . N 1, where ek ¼ 1= 2 if k ¼ 0, and ek ¼ 1 otherwise. The four point DCT matrix can be calculated from the above formula 0
0:5000 B 0:6533 A¼B @ 0:5000 0:2706
0:5000 0:2706 0:5000 0:6533
0:5000 0:2706 0:5000 0:6533
1 0:5000 0:6533 C C 0:5000 A 0:2706
(9.23)
The permutation matrix and TERMs factorized from matrix A are listed in Table 9.1.
0.0137 0.3066 1
0.6533 0.6533 0.5000 1
S 1 0 0 0.5307
1 0 0.8626
1 0.3933
9 Multiplierless Reversible Integer TDLT/KLT. . .
Table 9.1 Permutation matrix and TERMs factorized from 4-point DCT matrix P L U 0 1 0 0 1 1 0.2929 1 0 0 0 0.2346 1 1 0 0 0 1 0.4142 0.7654 1 0 0 1 0 0.2346 0 0.6934 1
1
195
196
J. Wu et al.
x0
y0 l10
u01 u02 u03
x1
u12 u13
x2 s30
s31 s32
l20
l21
P
y1 y2
l30
u23
l31
l32
x3
y3 S
U
L
Fig. 9.6 Four-point forward DCT implemented by multi-lifting
In this way, the matrix has been decomposed into TERMs, and diagonal elements of triangular matrices are equal to 1, so the 4-point DCT can be realized by multi-lifting as depicted in Fig. 9.6. From left to right, the input signals pass through S, U, L, and P in turn. If round the floating-point multiplication results to integers, then integer inputs will be transformed into integers. So in the inverse transform, original integers can be perfectly reconstructed as long as we subtract what was added as depicted in Fig. 9.5b.
3.2
Design of Reversible Integer to Integer TDLT/RKLT
In the proposed transform scheme, we use the modified matrix factorizing method, introduced above, to decompose the filtering matrix and DCT matrix into TERMs to realize reversible integer to integer transforms using the multi-lifting scheme. TDLT consists of pre- and post-filters with intact DCT between them. The preand post-filters are exact inverses of each other. The framework of TDLT [21] can be illustrated as Fig. 9.7. The general formula for the pre-filter [21] can be defined as: F¼
1 I 2 J
J I
I 0
0 V
I J
J ; I
(9.24)
where I and J are identity matrix and reversal identity matrix, respectively. Different types of TDLT can be derived with different matrices V. Two basic types of TDLT include time domain lapped orthogonal transform (TDLOT) and time domain lapped biorthogonal transform (TDLBT). The free-control matrix V is defined by the following two equations: VLOT ¼ JðCIIM=2 ÞT CIV M=2 J;
(9.25)
VLBT ¼ JðCIIM=2 ÞT DS CIV M=2 J;
(9.26)
9 Multiplierless Reversible Integer TDLT/KLT. . .
197
Fig. 9.7 Forward and inverse TDLT
a
Before pre-filtering
b
After pre-filtering Fig. 9.8 Pre-filtering effect
where CIIM=2 and CIV M=2 stand for the M/2 point type-II and type-IV DCT matrix respectively, DSp¼ffiffiffi diagfs; 1; ; 1gis a diagonal matrix where s is a scaling factor and we set s ¼ 2in our experiments. Now we take 2 2 filters to illustrate how they work. In this case, sub-matrices of F degenerate into matrices of one single 0 element with I ¼ J ¼ CIIM=2 ¼ CIV M=2 ¼ 1 and DS ¼ s. Letfxi gand fx i grepresent the input and output of the pre-filter, respectively, as depicted in Fig. 9.8. The 2 2 pre-filter operates on the two adjacent elements of neighboring blocks, so just x4 and x5 have been modified. The relationship between {x4 , x5 } and {x04 , x05 } can be obtained as follows. x0 4 x0 5
!
1 ¼ 2
1 1
1 1
1 0 0 s
1 1
1 1
x4 x5
!
(9.27)
198
J. Wu et al.
Fig. 9.9 Post-filtering effect
a
Before post-filtering
b
After post-filtering
1 x04 ¼ ½x4 þ x5 2
s ðx5
1 x05 ¼ ½x4 þ x5 þ s ðx5 2
x4 Þ ¼ x4 x4 Þ ¼ x5 þ
s
1 2
s
1 2
ðx5
x4 Þ
(9.28)
ðx5
x4 Þ
(9.29)
From (9.28) and (9.29) we can see that adjacent elements are modified by the pre-filter by enlarging their difference, aiming to lower the smaller-value pixel while increasing the larger-value pixel. Pre-filters with long-length basis not only enlarge the difference at boundaries from one block to another, but also make the pixels as homogenous as possible within each block. From another point of view, correlations between neighboring blocks have been reduced while correlations within one block have been increased. The post-filter is the exact inverse of the pre-filter, and n its o effect n on o the one_
dimensional vector has been depicted in Fig. 9.9, where x 0 i
_
and x i represent
the input and output of the post-filters. The post-filter decreases the difference between adjacent elements in neighboring blocks, aiming to increase the smallervalue pixel while lowering the larger-value pixel. For two-dimensional images, the post-filter aims to decrease blocking artifacts. The visual quality of post-filtered images will be shown in the experimental section, from which we will see that the post-filter performs very well. After pre-filtering, the input data will be transformed by DCT so that most energy can be concentrated in the low frequency region which will be beneficial to coding. Type-II DCT has been adopted in this paper. In the proposed method, we use the modified matrix factorization method [13] to decompose the filtering and DCT matrix into TERMs. This means that we do not need to search for other complex decomposition formations of the filtering matrix or DCT matrix, but simply factorize both of them to realize reversible integer-tointeger transform. TDLOT has largely reduced the blocking artifacts of DCT and greatly improved image compression but it cannot eliminate the blocking artifacts. TDLBT is
9 Multiplierless Reversible Integer TDLT/KLT. . .
199
constructed to modify TDLOT by inserting a diagonal matrix S in V acting as a scaling factor in order to further reduce the block artifacts. It should be pointed out that the determinant of the filter matrix F does not equal 1 and should be modified to satisfy det F ¼ 1 before factorizing TDLBT. In our method we normalize F in the following way: F ffiffiffiffiffiffiffiffiffiffiffiffiffi ; F ¼ p M det j Fj
(9.30)
where F is the original M M filter matrix, and F* is the normalized matrix. As soon as the filtering matrix and DCT matrix have been factorized into TERMs, RTDLT can be obtained, which includes reversible time domain lapped orthogonal transform (RTDLOT) and reversible time domain lapped biorthogonal transform (RTDLBT). If AF and ADCT denote the filtering matrix and DCT matrix, respectively, their factorizations will have the same format: AF ¼ PF LF UF SF
(9.31)
ADCT ¼ PDCT LDCT UDCT SDCT :
(9.32)
Therefore, we can obtain the multi-lifting structure of RTDLT if we replace the filters and DCT in Fig. 9.7 with the structure as in Fig. 9.6. PF , LF , UF , SF and PDCT , LDCT , UDCT , SDCT are TERMs decomposed from the matrices of Pre-filter and DCT, respectively. For RTDLT, because the PLUS factorizations only depend on the transform basis, we can calculate them in advance. Obviously, the best compression method should reduce redundancies in both the spatial and spectral dimensions for hyperspectral images. In our scheme, a spatial transform (RTDLT) is first applied. Then the spectral components of each spatial frequency band are decorrelated using a Reversible Integer Karhunen-Loeve Transform (RKLT).The dataset of hyperspectral images can be represented as below: X ¼ fX1 ; X2 ; X3 ; ; Xn gT ;
(9.33)
where the subscript n denotes the number of bands and Xn represent the sequence of the different spectral images. The W after RTDLT is represented as X
RTDLT
! W ¼ fW1 ;W2 ; W3 ; ; Wn gT
(9.34)
The covariance matrix Cw is defined as Cw ¼ E ¼
h W
m w ÞðW
X1 1 M ðW i M i¼0
mw ÞT
mw ÞðWi
i
m w ÞT
(9.35)
200
J. Wu et al.
m w ¼ EfW g ¼
X1 1 M Wi M i¼0
(9.36)
where M represents the number of the RTDLT coefficient vectors. Let li be the eigenvalues of matrixCw , with ei as the corresponding eigenvectors. When the eigenvalues li are arranged in descending order so that l1 l2 l3 li , the auto-correlations of the transformed signal vector are arranged in descending order. The transformation matrix TKLT is represented as TKLT ¼ fe1 ; e2 ; e3 ; ; en gT . The matrixY, a result of the KLT, is represented asY ¼ TKLT W. According to our proposal, the KLT matrix can be also realized by multiplierless multi-lift, based on matrix factorization. Integer approximation of the floating-point transform can be achieved if rounding operations are added, just as Eq. (9.2). Reversible integer KLT and TDLT are realized in the same way. The RKLT concentrates much of the energy into a single band, improving overall coding efficiency. KLT is the most efficient linear transform in the sense of energy compaction. The transform matrix can be obtained by calculating the eigenvectors of the covariance of the input data. To reduce the high computation complexity, lowcomplexity-KLT has been proposed by Penna et al. [32]. In our proposed method, integer reversible low-complexity-KLT (Low-RKLT) is designed for decorrelation in the spectral direction. The evaluation of the covariance matrix is simplified by sampling the input signal vectors. If using downsampling in KLT, then formula (9.36) should be rewritten as below: 0
1 MX1 mw ¼ 0 Wi : M i¼0 0
(9.37)
M denotes the number of the downsampled coefficients. In our experiments, 100:1 scaling is applied in the sampling process. As illustrated in [32], the performance of Low-RKLT is very similar to the fullcomplexity RKLT, but reduced the computation complexity significantly. The computational comparison will be given in Sect. 4. Reversible integer KLT and TDLT are realized in the same way. The integer transform can be realized by shifting and adding, only without any multiplier if the floating-point lifting coefficients are replaced by fractions, whose dominators are power of 2. For example, 15/64 ¼ 1/8 + 1/16 + 1/32 + 1/64 while 1/8, 1/16, 1/32 and 1/64 can be realized by shifting only. The multi-lifting coefficients of matrixces L, U, S decomposed from 4-point DCT have been tabulated in Table 9.2. Experimental results show that the multiplier-less DCT based on multi-lifting approximates the floating-point DCT very well. Further experiments to study the efficiency of multiplier-less RTDLT based on multi-lifting applied in lossy-to-lossless image compression will be discussed in the next section.
9 Multiplierless Reversible Integer TDLT/KLT. . . Table 9.2 Triangular TERM matrices with dyadic coefficients L U 1 1 75/256 1/64 167/256 15/64 1 1 39/128 167/256 53/128 49/64 1 1 1/2 15/64 0 89/128 1 1
201
S 1 0 1 0 0 1 17/32 221/256 101/256 1
4 Experimental Results and Discussion of RTDLT/RKLT-based Hyperspectral Image Compression A progressive hyperspectral image compression method is designed, based on RKLT and RTDLT, with higher compression ratio and better rate distortion (RD) performance compared with the 3D-DWT-based method. The flow graph is depicted in Fig. 9.10. In the spatial domain, 2-level RTDLT and 5-level DWT are adopted separately. In addition, as shown in Figs. 9.11 and 9.12, RTDLT transform coefficients are reorganized into a tree-structure [33] so as to be coded by waveletbased methods. In the coding method, 3DSPECK algorithm [34, 35] is applied, and JPEG2000-MC [36] has been applied in lossless. AVIRIS hyperspectral images [37], “Jasper Ridge” (scene 1), “Cuprite” (scene 3), “Lunar Lake” (scene 2), “Low Altitude” (scene 3) and “Moffett” (scene 1) (spatially cropped to 512 512, 224 bands) are used in our experiments to test the performance of different algorithms. The test images are coded with 16 bands in a group, so the entire image of 224 bands are divided into 14 groups. Figure 9.13 shows Jasper Ridge, Cuprite, Lunar Lake, Low Altitude and Moffett, respectively. Lossless compression performance is compared in Table 9.3, where two coding methods, 3DSPECK and JPEG2000-MC have been adopted, and transforms include asymmetric 3D-5/3DWT (3D-53DWT), 53DWT+RKLT and RTDLT+RKLT. Based on the same codec – 3DSPECK, we can see that our proposed RTDLT/RKLT performs 7.35 ~ 8.6% better than 3D-DWT and is comparable to 53DWT+RKLT. Lossy compression performance is given in Table 9.4. In our experiments, the performance of five types of transforms combined with 3D SPECK is compared at different bit rates, where 3D SPECK is carried with QccPack Version 0.56 [38]. It can be seen that our proposed RTDLT/RKLT consistently outperforms asymmetric 3D-97DWT by a large margin (up to 5 dB), combined with the same coding method. It also gives a gain of 0.38 ~ 0.69 dB compared with 53DWT+RKLT. Although the proposed method performs not as well as 97DWT+FloatKLT, it is capable of complete reversible transform while the latter cannot. Figure 9.14 depicts the rate distortion (RD) performance of different transform schemes using the same codec – 3D SPECK. From the above experimental results, we can conclude that the proposed compression method performs well in both lossy and lossless hyperspectral image compression. Among conventional transforms based on integer wavelet
202
J. Wu et al.
Fig. 9.10 Flow graph of the Proposed RTDLT/RKLT Compression Scheme
Fig. 9.11 Spatial coefficients reorganization for each band: (a) Original spatial coefficients distribution after RTDLT. (b) Reorganized coefficients as wavelet-liked subband structure
Fig. 9.12 Representations of block transform coefficients for Lena image: (a) Block representation after TDLT. (b) Wavelet-liked subband structure after coefficients reorganization
9 Multiplierless Reversible Integer TDLT/KLT. . .
203
Fig. 9.13 Original AVRIS images
(i.e. 5/3 DWT), the non-unitary transform gives the best lossless compression performance but will decrease the performance of lossy compression. Therefore, from a single lossless codestream, the decoder cannot obtain such good lossy results shown in Table 9.4. However, we don’t need to consider unitary problem by using
204
J. Wu et al.
Table 9.3 Lossless compression performance comparison (in bits per pixel per band, BPPPB) Codec 3D SPECK JPEG2000-MC Transform 3D-53DWT 53DWT + Low- RTDLT + low- 3D-53DWT 53DWT + LowRKLT RKLT RKLT Cuprite 5.32 4.97 4.95 5.44 5.07 Jasper 5.52 4.99 5.01 5.65 5.11 Lunar 5.29 4.97 4.96 5.42 5.08 Low 5.75 5.22 5.23 5.90 5.35 Moffett 6.67 6.11 6.11 6.86 6.32
Table 9.4 Lossy compression, SNR comparison (in dB) Transform methods/BPPPB 1 Cuprite Reversible 3D-53DWT 41.79 2D-53DWT + 1D-LowRKLT 43.77 2D-RTDLT + 1D-LowRKLT 44.61 Irreversible 3D-97DWT 43.24 2D-97DWT + 1D-LowFKLT 45.65 Jasper
Reversible
Irreversible Lunar
Reversible
Irreversible Low
Reversible
Irreversible Moffett
Reversible
Irreversible
0.75 40.64 43.07 43.91 41.93 44.65
0.5 38.78 42.21 43.03 39.89 43.54
0.25 34.72 39.29 40.13 35.77 40.46
0.125 30.73 34.33 35.16 31.62 35.38
3D-53DWT 2D-53DWT + 1D-LowRKLT 2D-RTDLT + 1D-LowRKLT 3D-97DWT 2D-97DWT + 1D-LowFKLT
34.80 38.63 39.45 36.27 40.54
33.01 37.57 38.43 34.19 39.28
30.35 35.75 36.52 31.39 37.05
25.52 30.7 31.45 26.49 31.77
21.06 25.14 25.93 22.07 26.06
3D-53DWT 2D-53DWT + 1D-LowRKLT 2D-RTDLT + 1D-LowRKLT 3D-97DWT 2D-97DWT + 1D-LowFKLT
42.75 44.55 45.34 44.34 46.42
41.68 43.85 44.61 43.04 45.44
39.89 43.05 43.75 41.07 44.36
35.84 40.42 41.01 36.98 41.51
31.63 35.61 36.10 32.65 36.67
3D-53DWT 2D-53DWT + 1D-LowRKLT 2D-RTDLT + 1D-LowRKLT 3D-97DWT 2D-97DWT + 1D-LowFKLT
34.36 38.04 38.68 35.65 39.49
32.82 37.06 37.66 33.87 38.34
30.51 35.65 36.19 31.46 36.72
26.31 31.84 32.22 27.21 32.64
22.43 26.89 27.33 23.28 27.61
3D-53DWT 2D-53DWT + 1D-LowRKLT 2D-RTDLT + 1D-LowRKLT 3D-97DWT 2D-97DWT + 1D-LowFKLT
41.24 43.36 43.99 41.83 44.22
39.61 41.69 42.36 40.18 42.62
37.04 39.27 39.83 37.64 40.03
31.55 33.49 34.15 32.57 34.21
26.28 27.64 28.34 27.32 28.37
uniform RTDLT framework to obtain embedded and high efficient lossy-to-lossless coding performance. In Table 9.5, we give a computation comparison between the low-complexity KLT and full-complexity KLT. From the table we can see that the low-complexity KLT has significantly reduced the computational time of full-complexity KLT.
9 Multiplierless Reversible Integer TDLT/KLT. . . Fig. 9.14 RD performance of different transform schemes combined with the same codec—3DSPECK
205
a 42 40 38 36
SNR (dB)
34 32 30 28 26
3D-53DWT 2D-53DWT+1D-LowRKLT 2D-RTDLT+1D-LowRKLT 3D-97DWT 2D-97DWT+1D-LowFKLT
24 22 20 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Bit Rate(bpppb)
Jasper Ridge
b
46 44 42
SNR (dB)
40 38 36 34
3D-53DWT 2D-53DWT+1D-LowRKLT 2D-RTDLT+1D-LowRKLT 3D-97DWT 2D-97DWT+1D-LowFKLT
32 30 28 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Bit Rate(bpppb)
Cuprite
Table 9.5 Computational time comparison (in seconds) Full-FKLT Low-FKLT Full-RKLT Low-RKLT Different KLTa Covariance computation 39.4375 0.2344 39.9844 0.2344 Spectral transform 71.8125 43.2813 98.2813 41.2031 a Full-FKLT full-complexity float KLT, Low-FKLT low-complexity float KLT, Full-RKLT fullcomplexity reversible KLT, Low-RKLT Low-complexity reversible KLT. The running times are measured on a Pentium IV PC at 3.20 GHz using C# Language
206
J. Wu et al.
Table 9.6 The performance comparison of various sampled KLT(SNR, in dB) Transform methods, sampling rates/BPPPB 1 0.75 0.5 Cuprite (512 512 224) Full-FKLT 43.81 43.06 42.22 Low-FKLT, 1/10 43.78 43.05 42.22 Low-FKLT, 1/100 43.76 43.02 42.19 Full-RKLT 43.83 43.13 42.26 Low-RKLT, 1/10 43.82 43.11 42.26 Low-RKLT, 1/100 43.77 43.07 42.21
0.25
0.125
39.37 39.38 39.33 39.34 39.35 39.29
34.43 34.42 34.37 34.35 34.36 34.33
Full-FKLT Low-FKLT, 1/10 Low-FKLT, 1/100 Full-RKLT Low-RKLT, 1/10 Low-RKLT, 1/100
Jasper (512 512 224) 38.67 37.57 35.78 38.67 37.57 35.77 38.68 37.57 35.79 38.56 37.48 35.66 38.61 37.55 35.71 38.63 37.57 35.75
30.75 30.73 30.71 30.68 30.69 30.70
25.22 25.20 25.19 25.05 25.06 25.14
Full-FKLT Low-FKLT, 1/10 Low-FKLT, 1/100 Full-RKLT Low-RKLT, 1/10 Low-RKLT, 1/100
Lunar (512 512 224) 44.53 43.81 43.01 44.53 43.81 43.02 44.53 43.81 43.02 44.55 43.84 43.05 44.55 43.84 43.05 44.55 43.85 43.05
40.49 40.49 40.48 40.44 40.46 40.42
35.74 35.73 35.69 35.62 35.63 35.61
However, there are also disadvantages, since the coding gain can be accompanied by lower speed. The complexity of RTDLT/RKLT is higher than that of 3D integer wavelet transforms. In Table 9.6, the performance comparison of various sampled 2D-53DWT+KLT is shown. The KLT compression performance degradation is negligible, when the scaling is 10:1~100:1 in the sampling process.
5 3D-RLT-based Hyperspectral Image Compression Following on from RTDLT/RKLT-based hyperspectral image compression, some readers might think intuitively whether RTDLT could be extended to apply to a spectral domain. Certainly, the RTDLT can be extended to 3D reversible integer lapped transform (3D-RLT), which is realized by cascading three 1D RLT along spatial and spectral dimensions. RLT is derived from reversible integer time-domain lapped transform (RTDLT) which was introduced in Sect. 3. Compared with RTDLT/RKLT, the new 3D block transform – 3D-RLT is much less complex. We have extended RLT to three dimensions, and designed a 3D-RLT-based compression algorithm for hyperspectral images as follows.
9 Multiplierless Reversible Integer TDLT/KLT. . .
207
Output
Reversible Lapped Transform (RLT) on the spectral dimension
Reversible Integer DCT (RDCT)
Reversible Integer DCT (RDCT)
Reversible Integer Pre-Filter (R-P-F)
Reversible Integer DCT (RDCT)
Reversible Integer Pre-Filter (R-P-F)
Reversible Integer DCT (RDCT)
Reversible Integer Pre-Filter (R-P-F)
Input
Band n
Band 3
Band 2
Band 1
…………
Fig. 9.15 The diagram of reversible lapped transform (RLT) on the spectral dimension
It is necessary to perform 2D RLT on each band of the entire image along the horizontal and vertical dimensions, followed by 1D RLT along the spectral dimension, as shown in Fig. 9.15. RLT is realized by performing a reversible integer prefilter (R-P-F) and reversible integer DCT (RDCT) successively. In the spatial domain, RLT is performed on the unit of a block. Here, we adopt blocks with size 8 8. The test images are coded with 16 bands in a group in the spectral dimension. R-P-F can be turned on or off depending on the performance requirement and complexity limitation. If R-P-F is turned off, RLT converts to RDCT.
6 Experimental Results of 3D-RLT-based Hyperspectral Image Compression To demonstrate the validity of the proposed 3D-RLT coding scheme, we again conducted experiments on the AVIRIS hyperspectral images [37]. 3D floating-point 9/7-tap biorthogonal WT (3D-97WT), 3D integer 5/3-tap WT (3D-53WT), 3D floating-point DCT (3D-FDCT), 3D-RDCT, 3D floating-point TDLT (3D-FLT) and 3D-RLT have been tested. It has been proved that the asymmetric 3D-WT (anisotropic wavelet) performs better than the symmetric 3D-WT (isotropic wavelet) [39]. In our experiment, we adopt asymmetric 3D-WT. RLT coefficients are approximated by hardware-friendly dyadic values which realize operations in the transforming process be realized by only shifts and additions. All the transform methods are combined with the same coding method, namely the 3D set partitioned embedded block (3D-SPECK) codec [35], to ensure a fair comparison. Lossy compression results, based on the criterion of signal-to-noise ratio (SNR), are given in Table 9.7, from which we can see that at most bit-rates 3D-RLT performs better than all the other reversible transform methods except the floating-point transforms. Although the floating-point transforms produce the better results, 3D-RLT can be applied in progressive lossy-to-lossless compression.
208
J. Wu et al.
Table 9.7 Lossy compression performance (SNR, in dB) Transform methods/BPPPB 1 Cuprite Reversible 3D-53DWT 41.79 3D-RDCT 42.86 3D-RLT 42.76 2D-RLT + 1D-RDCT 42.76 Irreversible 3D-97DWT 43.24 2D-97DWT + 1D-FDCT 43.48 2D-97DWT + 1D-FLT 43.56 3D-FDCT 43.34 3D-FLT 43.53 Jasper
Reversible
Irreversible
Lunar
Reversible
Irreversible
Low
Reversible
Irreversible
Moffett
Reversible
Irreversible
0.75 40.64 41.65 41.81 41.66 41.93 42.23 42.39 42.02 42.35
0.5 38.78 39.77 40.25 39.92 39.89 40.24 40.64 40.01 40.67
0.25 34.72 35.57 36.71 35.87 35.77 36.11 36.92 35.73 36.98
0.125 30.73 31.49 32.59 31.89 31.62 31.95 32.68 31.59 32.83
3D-53DWT 3D-RDCT 3D-RLT 2D-RLT + 1D-RDCT 3D-97DWT 2D-97DWT + 1D-FDCT 2D-97DWT + 1D-FLT 3D-FDCT 3D-FLT
34.80 35.78 35.73 35.79 36.27 36.32 36.35 36.09 36.24
33.01 33.92 34.07 34.03 34.19 34.41 34.52 34.13 34.49
30.35 31.29 31.80 31.47 31.39 31.64 32.00 31.39 32.02
25.52 26.48 27.39 26.79 26.49 26.81 27.45 26.53 27.43
21.06 22.10 23.09 22.51 22.07 22.43 23.01 22.06 23.03
3D-53DWT 3D-RDCT 3D-RLT 2D-RLT + 1D-RDCT 3D-97DWT 2D-97DWT + 1D-FDCT 2D-97DWT + 1D-FLT 3D-FDCT 3D-FLT
42.75 43.77 43.66 43.68 44.34 44.44 44.51 44.26 44.45
41.68 42.59 42.72 42.62 43.04 43.20 43.36 42.96 43.31
39.89 40.65 41.14 40.85 41.07 41.31 41.66 40.91 41.59
35.84 36.44 37.59 36.87 36.98 37.28 38.01 36.69 37.92
31.63 32.16 33.38 32.67 32.65 33.01 33.71 32.33 33.68
3D-53DWT 3D-RDCT 3D-RLT 3D-97DWT 2D-97DWT + 1D-FDCT 2D-97DWT + 1D-FLT 3D-FDCT 3D-FLT
34.36 35.04 35.07 35.65 35.55 35.66 35.31 35.51
32.82 33.48 33.62 33.87 33.92 34.16 33.65 33.92
30.51 31.13 31.65 31.46 31.62 32.03 31.22 31.76
26.31 26.92 27.72 27.21 27.46 28.04 26.94 27.79
22.43 23.01 23.97 23.28 23.56 24.23 22.95 23.91
3D-53DWT 3D-RDCT 3D-RLT 3D-97DWT 2D-97DWT + 1D-FDCT 2D-97DWT + 1D-FLT 3D-FDCT 3D-FLT
41.24 40.72 37.82 41.83 40.91 37.96 40.77 37.83
39.61 39.23 36.69 40.18 39.46 36.85 39.26 36.75
37.04 36.91 34.88 37.64 37.26 35.13 36.96 34.98
31.55 31.97 31.06 32.57 32.41 31.37 32.03 31.24
26.28 26.73 26.59 27.32 27.18 26.93 26.71 26.79
9 Multiplierless Reversible Integer TDLT/KLT. . .
209
Fig. 9.16 Spectral profiles of various images
3D-RLT performs better than 3D-53DWT by about 0.71~2.03 dB, and better than 3D-97DWT in most cases. However, for the Moffett image, it is clear in comparison with 3D-DWT the performance of 3D-RLT is degraded significantly. This is because Moffett’s spectral profile is much smoother than that of other images, as shown in Fig. 9.16. It is well-known, either DCT or lapped transform can be considered to be a kind of multi-band filters, so DCT and lapped transform are not suitable in the spectral dimension of this kind of images which has simple frequency characteristic. The lossless compression performance of reversible integer transforms combined with 3DSPECK is compared in Table 9.8 based on the criterion of bit per pixel per band (bpppb). From Table 9.8 we can see that 3D-RLT is competitive with 3D-53DWT and better than 3D-RDCT in most cases.
210 Table 9.8 Lossless compression performance (in bits per pixel per band, BPPPB) Images 3D-53DWT 3D-RDCT Cuprite 5.32 5.38 Jasper 5.52 5.65 Lunar 5.29 5.36 Low 5.75 5.87 Moffett 6.67 7.31
J. Wu et al.
3D-RLT 5.33 5.63 5.32 5.83 7.61
7 Conclusions In this chapter, we first presented a new transform scheme of RTDLT/RKLT for the compression of hyperspectral images with a performance that is competitive or even superior to that of state-of-the-art conventional transform-based techniques. The proposed RTDLT/RKLT scheme can realize reversible integer-to-integer transform, and so it can be applied in progressive lossy-to-lossless compression combined with zero-block-based bit-plane coders. The new transforming scheme has some advantages such as multiplierless computing, flexibility in processing, and parallel implementation and so on. In addition, it also preserves the desirable features of the lifting scheme, such as low memory request, in-place calculation, and perfect reconstruction. While high performance means high complexity, compared with 1D wavelet transform, RKLT still has much higher complexity in the spectral dimension. We also presented a simple but efficient compression algorithm based on 3DRLT for hyperspectral images. In other words, the 1D reversible integer lapped transform is applied on all dimensions of hypespectral images. Our experiments showed that 3D-RLT can defeat the integer 5/3 wavelet easily in a lossy application. At the same time, the lossless performance of 3D-RLT is near that of 5/3 WT. In addition, the 3D-RLT-based method can be simplified if much lower complexity is required. When the pre/post filters are cancelled, the RLT can be converted into RDCT. At this point, it should be noted that in our experiments 5/3 WT is unitary in lossy mode in order to obtain better performance, so the experimental results of 5/3 WT do not support completely progressive transforms. If, in order to support progressive transforms completely and, at the same time, to ensure the best lossless performance, 5/3 WT should be non-unitary, although that will degenerate the lossy performance. In contrast, the unitary wavelet transform will cause the degeneration of lossless compression. However, our proposed methods do not suffer this problem, and both RTDLT/RKLT and 3D-RLT supports progressive coding completely. No transform provides a perfect result and is only one of many topics relating to image compression. Although our proposed RTDLT/RKLT and 3D-RLT methods are suitable for the lossy-to-lossless compression of hyperspectral images (actually, RTDLT is also suitable for the compression of natural images [40]), they are not suitable to transform 3D medical images. The reason is because DCT-based lapped
9 Multiplierless Reversible Integer TDLT/KLT. . .
211
transform is not suitable for medical images with a smooth texture. In this case, the wavelet transform is rather good. In addition, if a given application is only for 3D lossless compression of hyperspectral images, it would be better to use the lookup table (LUT) coding method [41] which has low complexity and high lossless performance.
References 1. J. Xu, Z. Xiong, S. Li, and Y. Zhang, “3-D embedded subband coding with optimal truncation (3-D ESCOT),” Applied and Computational Harmonic Analysis, vol.10, pp.290–315, May 2001. 2. J. E. Fowler and D. N. Fox, “Embedded wavelet-based coding of three dimensional oceanographic images with land masses,” IEEE Transactions on Geoscience and Remote Sensing, vol.39, no.2, pp.284–290, February 2001. 3. X. Tang, W. A. Pearlman, and J. W. Modestino, “Hyperspectral image compression using three-dimensional wavelet coding,” in proceedings SPIE, vol.5022, pp.1037–1047, 2003. 4. Q. Du and J. E. Fowler, “Hyperspectral Image Compression Using JPEG2000 and Principal Component Analysis,” IEEE Geoscience and Remote sensing letters, vol.4, pp.201–205, April 2007. 5. P. L. Dragotti, G. Poggi, and A. R. P. Ragozini, “Compression of multispectral images by three-dimensional SPIHT algorithm,” IEEE Geoscience and Remote sensing letters, vol.38, no. 1, pp. 416–428, January 2000. 6. B. Penna, T. Tillo, E. Magli, and G. Olmo, “Transform Coding Techniques for Lossy Hyperspectral Data Compression,” IEEE Geoscience and Remote sensing letters, vol.45, no.5, pp.1408–1421, May 2007. 7. B. Penna, T. Tillo, E. Magli, and G. Olmo, “Progressive 3-D coding of hyperspectral images based on JPEG 2000,” IEEE Geoscience and Remote sensing letters, vol.3, no.1, pp.125–129, January 2006. 8. W. Sweldens, “The lifting scheme: A construction of second generation wavelet,” in SIAM Journal on Mathematical Analysis, vol.29, pp.511–546, 1997. 9. I. Daubechies and W. Sweldens, “Factoring wavelet transforms into lifting steps,” The Journal of Fourier Analysis and Applications, vol.4, pp.247–269, 1998. 10. A.Bilgin, G. Zweig, and M. W. Marcellin, “Three-dimensional image compression using integer wavelet transforms,” Applied Optics, vol.39, pp.1799–1814, April 2000. 11. Z. Xiong, X. Wu, S. Cheng, and J. Hua, “Lossy-to-Lossless Compression of Medical Volumetric Data Using Three-Dimensional Integer Wavelet Transforms,” IEEE Transactions on Medical Imaging, vol.22, no.3, pp.459–470, March 2003. 12. P. Hao and Q. Shi, “Reversible integer KLT for progressive-to-lossless compression of multiple component images,” in Proceedings IEEE International Conference Image Processing (ICIP’03), Barcelona, Spain, pp.I-633–I-636, 2003. 13. L. Galli and S. Salzo, “Lossless hyperspectral compression using KLT,” IEEE International Geoscience and Remote Sensing Symposium, (IGARSS2004), vol.1, pp.313–316, September 2004. 14. C. Kwan, B. Li, R. Xu, X. Li, T. Tran, and T. Nguyen, “A Complete Image Compression Method Based on Overlapped Block Transform with Post-Processing,” EURASIP Journal on Applied Signal Processing, pp.1–15, January 2006. 15. P. List, A. Joch, J. Lainema, G. Bjontegaard, M. Karczewicz, “Adaptive deblocking filter,” IEEE Transactions on Circuits and Systems for Video Technology, vol.13, no.7, pp.614–619, July 2003.
212
J. Wu et al.
16. Xiong ZX, Orchard MT, Zhang YQ, “A deblocking algorithm for JPEG compressed images using overcomplete wavelet representations,” IEEE Transactions on Circuits and Systems for Video Technology, vol.7, no.2, pp.433–437, April 1997. 17. P. Cassereau, “A New Class of Optimal Unitary Transforms for Image Processing”, Master’s Thesis, Massachusetts Institute of Technology, Cambridge, MA, May 1985. 18. H. S. Malvar, “Lapped transforms for efficient transform/subband coding”, IEEE Transactions on Acoustics, Speech, and Signal Processing, pp.969–978, ASSP-38. 1990. 19. C.W. Lee and H. Ko, “Arbitrary resizing of images in DCT domain using lapped transforms”, Electronics Letters, vol.41, pp.1319–1320, November 2005. 20. T. D. Tran, J. Liang, and C. Tu, “Lapped transform via time-domain pre- and post-processing,” IEEE Transactions on Signal Processing, vol.51, no.6, pp.1557–1571, January 2003. 21. Chengjie Tu, and Trac D. Tran, “Context-based entropy coding of block transform coefficients for image compression,” IEEE Transactions on Image Processing, vol.11, no.11, pp. 1271–1283, January 2002. 22. Jie Liang, Trac D. Tran, “Fast Multiplierless Approximations of the DCT with the Lifting Scheme,” IEEE Transactions on Signal Processing, vol.49, no.12, pp.3032–3044, December 2001. 23. http://www.microsoft.com/whdc/xps/hdphotodpk.mspx 24. S. Srinivasan, C. Tu, S. L. Regunathan, and G. J. Sullivan, “HD Photo: a new image coding technology for digital photography,” in Proceedings SPIE Applications of Digital Image Processing XXX, San Diego, vol.6696, pp.66960A, August 2007. 25. C. Tu, S. Srinivasan, G. J. Sullivan, S. Regunathan, and H. S. Malvar, “Low-complexity hierarchical lapped transform for lossy-to-lossless image coding in JPEG XR//HD Photo,” in Proceedings SPIE Applications of Digital Image Processing XXXI, San Diego,vol.7073, pp.70730 C1-12,August 2008. 26. Y. Chen, S. Oraintara, and T. Nguyen, “Integer discrete cosine transform (IntDCT),” in Proceedings 2nd International Conference Information and Communication of Signal Processing, December 1999. 27. G.C.K. Abhayaratne, “Reversible integer-to-integer mapping of N-point orthonormal block transforms,” Signal Processing, vol.87, no.5, pp.950–969, 2007. 28. J. Li, “Reversible FFT and MDCT via matrix lifting,” in Proceedings IEEE International Conference on Acoustics, Speech, and Signal Processing, vol. 4, pp. iv-173–iv-176, May 2004. 29. L.Z. Cheng, G.J. Zhong, and J.S. Luo, “New family of lapped biorthogonal transform via lifting steps,” IEE Proceedings -Vision, Image and Signal Processing, Vol. 149, no. 2, pp. 91–96, April 2002. 30. P. Hao and Q. Shi, “Matrix Factorizations for Reversible Integer Mapping,” IEEE Transactions on Signal Processing, vol.42, no.10, pp. 2314–2324,October 2001. 31. K. R. Rao and P. Yip, Discrete Cosine Transform: Algorithms, Advantages, Applications. New York: Academic, 1990. 32. B. Penna, T. Tillo, E. Magli, and G. Olmo, “Transform Coding Techniques for Lossy Hyperspectral Data Compression,” IEEE Geosciense and Remote Sensing, vol. 45, no. 5, pp. 1408–1421, May 2007. 33. Z. Xiong, O. Guleryuz, and M. T. Orchard, “A DCT-based embedded image coder”. IEEE Signal Processing Letters, vol. 3, pp. 289–290, 1996. 34. X. Tang and W. A. Pearlman, “Three-Dimensional Wavelet-Based Compression of Hyperspectral Images,” in Hyperspectral Data Compression, pp.273–308, 2006. 35. Jiaji Wu, Zhensen Wu, Chengke Wu, “Lossly to lossless compression of hyperspectral images using three-dimensional set partitioning algorithm,” Optical Engineering, vol.45, no.2, pp.1–8, February 2006. 36. Information Technology—JPEG 2000 Image Coding System—Part 2: Extensions, ISO/IEC 15444–2, 2004. 37. http://aviris.jpl.nasa.gov/html/aviris.freedata.html 38. http://qccpack.sourceforge.net/
9 Multiplierless Reversible Integer TDLT/KLT. . .
213
39. E. Christophe, C. Mailhes and P. Duhamel, “Hyperspectral image compression: adapting SPIHT and EZW to Anisotropic 3-D Wavelet Coding”, IEEE Transactions on Image Processing, vol.17, no.12, pp.2334–2346, 2008. 40. Lei Wang, Jiaji Wu, Licheng Jiao, Li Zhang and Guangming Shi, “Lossy to Lossless Image Compression Based on Reversible Integer DCT. IEEE International Conference on Image Processing 2008 (ICIP2008), pp.1037–1040, 2008. 41. J. Mielikainen, “Lossless compression of hyperspectral images using lookup tables,” IEEE Signal Processing Letters, vol.13, no.3, pp.157–160, 2006.
Chapter 10
Divide-and-Conquer Decorrelation for Hyperspectral Data Compression Ian Blanes, Joan Serra-Sagrista`, and Peter Schelkens
Abstract Recent advances in the development of modern satellite sensors have increased the need for image coding, because of the huge volume of such collected data. It is well-known that the Karhunen-Loeˆve transform provides the best spectral decorrelation. However, it entails some drawbacks like high computational cost, high memory requirements, its lack of component scalability, and its difficult practical implementation. In this contributed chapter we revise some of the recent proposals that have been published to mitigate some of these drawbacks, in particular, those proposals based on a divide-and-conquer decorrelation strategy. In addition, we provide a comparison among the coding performance, the computational cost, and the component scalability of these different strategies, for lossy, for progressive lossy-to-lossless, and for lossless remote-sensing image coding.
1 Introduction When coding a hyperspectral image in a lossy or progressive lossy-to-lossless way, it is common to employ a spectral decorrelating transform followed by a traditional transform coder. In this regard, the Karhunen-Loeˆve Transform (KLT) and its derivatives are the transforms that provide some of the best results [20]. Yet, the KLT has a very high computational cost that hinders its adoption in many situations.
I. Blanes (*) • J. Serra-Sagrista` Universitat Auto`noma de Barcelona, E-08290 Cerdanyola del Valle`s (Barcelona), Spain e-mail: [emailprotected]; [emailprotected] P. Schelkens Department of Electronics and Informatics, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium Interdisciplinary Institute for Broadband Technology, Gaston Crommenlaan 8, b102, B-9050 Ghent, Belgium e-mail: [emailprotected] B. Huang (ed.), Satellite Data Compression, DOI 10.1007/978-1-4614-1183-3_10, # Springer Science+Business Media, LLC 2011
215
216
I. Blanes et al.
c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15
Fig. 10.1 Example of divide-and-conquer strategy for 15 spectral components
Recently, a variety of divide-and-conquer strategies have been proposed to ease the computational requirements of the KLT. These strategies rely on the fact that the KLT has a quadratic cost O(n2), but that if the full transform is approximated by a set of smaller ones, then the computational cost becomes a fraction of the original cost. The rationale behind these strategies is that only spectral components with high covariances are worth decorrelating, since in the other cases, coding gains are negligible. Hence, in these strategies, a transform is decomposed into a set of smaller transforms that provide decorrelation where it is needed and the other regions are neglected. One of such transform decompositions is shown in Fig. 10.1. In this example, first, local decorrelation is provided by three KLT clusters, and then, some outputs of each clusters are processed again together to achieve global decorrelation. Note that due to the properties of the KLT, output components are by convention arranged in descending order according to their variance, therefore, components with low variances are discarded by not selecting the last outputs from each cluster. This chapter reviews divide-and-conquer strategies for the KLT, and provides some insights into the details found in their practical implementation. Additionally, it includes an experimental evaluation and comparison of the various strategies in the context of the compression of hyperspectral remote-sensing images. We note that some other recent approaches to alleviate some of the issues of the KLT are not addressed here. In particular, the reader should be aware that coding gains can be improved by not assuming Gaussian sources, and finding the Optimal Spectral Transform with an Independent Component Analysis (ICA)-based algorithm [3, 4], which could also be pre-trained [2, 5]; or that, to overcome the high computational cost of the KLT, the Discrete Cosine Transform (DCT) was proposed [1] – it assumes a Toeplitz matrix as data covariance matrix. However, the DCT has a poor performance as spectral decorrelator [20]. Similar approaches to reduce the computational cost are the fast Approximate Karhunen-Loeˆve Transform (AKLT) [19] and the AKLT2 [21], which extend the DCT with first and second order perturbations. The rest of this chapter is organized as follows: Sect. 10.2 provides an overview of the Karhunen-Loeˆve Transform. In Sect. 10.3 we present the different variants of divide-and-conquer strategies for spectral decorrelation. Section 10.4 provides some experimental results for compression of satellite data using the different decorrelation strategies. Finally, Sect. 10.5 contains some conclusions.
10
Divide-and-Conquer Decorrelation for Hyperspectral Data Compression
217
2 The Karhunen-Loeˆve Transform This section describes the KLT and its practical application when used as part of a divide-and-conquer strategy. For an image of N spectral components, each one with a mean required to be zero, the KLT is defined as follows: let X be a matrix that has N rows, one for each source, and M columns, one for each spatial location. Then, Y , the outcome after applying the transform, is computed as Y ¼ KLTSX ðXÞ ¼ QT X;
(10.1)
where SX ¼ ð1=MÞXXT is the covariance matrix of X, and Q is the orthogonal matrix obtained from the Eigenvalue Decomposition (ED) of SX, i.e., SX ¼ QLQ 1 , with L ¼ diag(l1, . . ., lN), jl1j jl2j . . . jlNj . Being SX an Hermitian matrix, according to the spectral theorem, such a decomposition always exists. The covariance matrix of Y is the diagonal matrix L (i.e., SY ¼ ð1=MÞYY T ¼ ð1=MÞQT XXT Q ¼ QT SX Q ¼ L), and l1, . . ., lN are the variances of each component after the transform. The ED of the covariance matrix is usually performed using an iterative procedure that converges to the solution, such as the algorithm based on a preprocessing tridiagonalization by Householder transformations, followed by iterative QR decompositions [13]. Other diagonalization procedures with lower complexities exist [10], at the expense of higher implementation difficulties. It is highly recommended to rely on one of the existing libraries for this purpose, as several numerical stability issues appear in the process that may lead to a non-converging algorithm. Note that this transform is dependent on SX and hence different for each input. For this reason, the specific transform matrix QT used has to be signaled as side information, so that the decoder is able to invert the transformation.
2.1
Centering and Covariance Particularities
As noted in its definition, the KLT has to be applied on image components of zero mean. It is rare to find image components with zero mean, thus, a variable change is usually applied to address this issue, i.e., 1 0 M P 1 xð1; jÞ C BM C B j¼1 C B .. Cð1; . . . ; 1Þ: (10.2) X0 ¼ X B C B . C B P M A @1 xðN; jÞ M j¼1
218
I. Blanes et al.
As the KLT does not affect the mean of components, this change of variable is only required once per image, regardless of whether a sequence of transforms or just one transform is applied. Some particularities of the covariance calculation are also worth noting. With large enough spatial sizes, image subsampling in the KLT has been known to substantially reduce the training cost of the KLT [20]. The aim is to use only a small sample of the image spatial locations to compute the covariance matrix. A subsampling factor of r ¼ 0. 01 (1%) is known to provide almost the same compression performance results, and effectively reduce transforms training cost, leaving the application of the transform matrix – the QTX matrix multiplication – as the main source of computational cost. To select sampling locations, some sort of Pseudo-Random Number Generator (PRNG) is required. High quality PRNG have high cost, but, in this case, with a poor generator, results stay similar. A very fast Park-Miller PRNG (10.3) can be used, only taking four operations for each random number. Y n ¼ M Xn div ð232
5Þ
Xn ¼ ð279470273 Xn 1 Þ mod ð232
5Þ
(10.3)
Another particularity of the covariance matrix calculation appears when the KLT is used in a divide-and-conquer strategy, where one or more components are transformed more than once. When a component has been previously transformed and has to be transformed again, its variance does not have to be computed, as it is already available as one of the eigenvalues li from the previous transform.
2.2
Spatial Subdivision
Spatial subdivision – also known as segmentation or blocking – might be used in this context for several reasons, e.g., to ease memory requirements or to limit the impact of data loss. Some issues are here discussed in relation to its use. If a transform is applied in multiple spatial blocks, or an image spatial size is small enough, then the computational cost for the ED step stops being negligible and might mean an important fraction of the total computational cost. To address this issue, improved methods of ED can be used [10, 14]. It is also worth noting that the size of the transform side information is only dependent on the spectral size, thus with multiple spatial divisions or small spatial sizes, the size of the transform side information becomes more relevant, as the total amount of data per transform becomes smaller. Covariance subsampling is also affected by spatial subdivision for the same reasons; as the spatial size decreases, the number of sampled pixels also decreases. Hence, subsampling factors have to be increased –i.e., a larger number of samples has to be used–, to still produce good estimations. When spatial sizes are very small, covariance subsampling might even
10
Divide-and-Conquer Decorrelation for Hyperspectral Data Compression
219
be more costly than a regular covariance calculation because of the small overhead of sample selection becoming noticeable. Finally, in case of multiple spatial blocks, blocking artifacts might also appear on the spatial boundaries and have to be taken into consideration when evaluating specific transforms.
2.3
The Reversible KLT
Given an almost exactly-reversible floating-point KLT matrix, like the one obtained in the previous steps, one can adapt it to be used in a lossless way. This adaptation is performed factorizing the transform matrix into an equivalent sequence of Elementary Reversible Matrices (ERMs) [9]. Each of these ERMs can later be applied in an approximate way with respect to the original multiplication, which maps integers to integers, and is fully reversible. In addition to the reversibility, as the transformed coefficients are integers instead of floating-point values, the number of bitplanes to be encoded later is reduced. A factorization known as the Reversible KLT (RKLT) with quasi-complete pivoting [12] is described now. The factorization is based on the one introduced in [16], and improves it by minimizing the adaptation differences with the original lossy transform. The factorization is as follows: given an orthogonal KLT matrix QT, an iterative process over Ai is applied, starting with A1 ¼ QT. 1. A permutation matrix Pi is selected such that (PiAi)(i, N) 6¼ 0. 2. The matrix Si is computed as Si ¼ I si ek ei T , where em is the m-th vector of the standard basis, and si ¼
ðPi Ai Þðj;iÞ 1 ðPi Ai Þðj;kÞ .
Indices j and k are selected so that si is
minimal over i j N, i þ 1 k N. 3. Gaussian elimination is applied to PiAiSi, obtaining the Gaussian elimination matrix Li that guarantees (LiPiAiSi)(k, i) ¼ 0 for all k > i. 4. Ai þ 1 is set to the product LiPiAiSi, and a new iteration is started. After N 1 iterations, AN is an upper triangular matrix, from now on called U, where all diagonal elements are 1 except the last one which might be 1. Then all the partial results obtained in the previous step are merged, S
1
¼
N Y1
Sk ;
(10.4)
T
(10.5)
k¼1
L
1
¼ LN
1
ðPN 1 LN 2 PN
1
Þ ðP2 L1 P2 T Þ;
220
I. Blanes et al.
PT ¼
N Y1
PN k ;
(10.6)
k¼1
and finally the factorization is reached: A ¼ PLUS;
(10.7)
where P is a permutation matrix, and L and S are lower triangular matrices with ones in the diagonal. Once the matrix is factorized, the transform is applied in the forward and inverse directions using interleaved rounding operations within each specially crafted matrix multiplication. The rounding operator —denoted by []— is defined as the function that rounds its argument to the closest integer value, with half-way values rounded to the nearest even value. For lower triangular matrices, the multiplication – for instance for Y ¼ LX– is as follows: 1 1 0 x1 y1 . C B .. C B hP .. i C B . C B C B B C B x þ C C¼B i 1j