**The Wavelet Transform as a Classification Criterion Applied to Improve Compression of Hyperspectral Images**

Daniel Acevedo and Ana Ruedin *Departamento de Computación Facultad de Ciencias Exactas y Naturales Universidad de Buenos Aires Argentina*

### **1. Introduction**

526 Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology

combined transforms could compensate for the drawbacks of each other, resulting in

[1] S.P. Mohanty, et al., "A Dual Watermarking Technique for Images", Proc.7th ACM

[2] Yusnita Yusof and Othman O. Khalifa , "Digital Watermarking For Digital Images Using

[3] Potdar, V., S. Han and E. Chang, 2005. "A Survey of Digital Image Watermarking

[4] Chan, C. and L. Cheng, 2004. "Hiding Data in Images by Simple LSB Substitution",

[5] Wang, R., C. Lin and J. Lin, " Copyright protection of digital images by means of

[6] G. Langelaar, I. Setyawan, R.L. Lagendijk, "Watermarking Digital Image and Video Data", in IEEE Signal Processing Magazine, Vol 17, pp 20-43, September 2000. [7] Ingemar J.Cox, Matt L. Miller and Jeffrey A. Bloom, "Watermarking Applications and their properties", Int. Conf. On Information Technology'2000, Las Vegas, 2000. [8] Gerhard C. Langelaar, Iwan Setyawan, and Reginald L. Lagendijk, "Watermarking

[9] Maurice Mass, Ton Kalker, Jean-Paul M.G Linnartz, Joop Talstra, Geert F. G. Depovere,

[10] Fabien A.P. Petitcolas, " Watermarking Schemes Evaulation", IEEE Signal Processing

[11] Technical Report, submitted to The Scientific and Technical Research Council of Turkey

[12] Jean François Delaigle, " Protection of Intellectual Property of Images by perceptual

[13] Ingemar J. Cox, Joe Kilian, Tom Leighton, and Talal Shamoon, "Secure Spread

[14] Mitchell D. Swanson, Mei Kobayashi, and Ahmed H. Tewfik, "Multimedia Data-

[15] Mitchell D. Swanson, Bin Zhu, and Ahmed H. Tewfik, "Transparent Robust Image Watermarking", 1996 SPIE Conf. on Visual Communications and Image Proc. [16] Christine I. Podilchuk and Wenjun Zeng, " Image-Adaptive Watermarking Using Visual

[17] Raymond B. Wolfgang, Christine I. Podilchuk and Edward J. Delp, "Perceptual

[18] Sergio D. Servetto, Christine I. Podilchuk, Kannan Ramchandran, "Capacity Issues in

Conference on Image Processing (ICIP), Chicago, IL, October 1998.

Wavelet Transform", Proc 2007 IEEE conference, pp 665-669

Data/Image Coding, Compression, and Encryption, USA.

Informatics, pp: 709-716, Perth, Australia.

Signal Processing Magazine, September 2000.

(Tübitak) under project EEEAG 101E007, April 2002

Sciences, Universite Catholique de Louvain, Belgique.

Magazine, September 2000

1673-1687, (1997).

No. 6, June 1998.

Pattern Recognition, 37(3):469-474.

International Multimedia Conference, ACM-MM'99, Part 2, pp 49-51, Orlando,

Techniques", in Proc. of the IEEE International Conference on Industrial

frequency domain watermarking," Proc. of the SPIE Conference On Mathematics of

Digital Image and Video Data", IEEE Signal Processing Magazine, September 2000.

and Jaap Haitsma, " Digital Watermarking for DVD Video Copy Protection", IEEE

Watermarking", Ph.D Thesis submitted for the degree of Doctor of Applied

Spectrum Watermarking for Multimedia", IEEE Trans. on Image Processing, 6, 12,

Embedding and Watermarking Technologies", Proceedings of the IEEE., Vol. 86,

Models", IEEE Journal of Selected Areas in Communications, Vol.16, No.4, May 1998.

Watermarks for Image and Video", Proceedings of the IEEE, Vol. 87, No. 7, July 1998.

Digital Image Watermarking", In the Proceedings of the IEEE International

effective watermarking.

USA, Oct. 1999

**5. References** 

Satellites continually feeding images to their base, pose a challenge as to the design of compression techniques to store these huge data volumes. We aim at lossless compression of hyperspectral images having around 200 bands, such as AVIRIS images. These images consist of several images (or bands) obtained by filtering radiation from the earth at different wavelengths. Compression is generally achieved through reduction of spatial as well as spectral correlations.

Most of the hyperspectral compressors are prediction–based. Since spectral correlation is usually high (much more higher than spatial correlation) pixels are predicted with other pixels in an adjacent band (rather than other pixels surrounding the one to be predicted). SLSQ (Rizzo et al., 2005), a low-complexity method designed for hyperspectral image compression, performs a simple prediction for each pixel, by taking a constant times the same pixel in the previous band. The constant is calculated by least squares over 4 previously encoded neighboring pixels. SLSQ–OPT version of SLSQ performs one band look–ahead to determine if the whole band is better compressed this way or with intraband prediction, while in the SLSQ–HEU version this decision is taken by an offline heuristic. CCAP (Wang et al., 2005) predicts a pixel with the conditional expected value of a pixel given the context. The expected value is calculated over coded pixels having matching (highly correlated) contexts. Slyz and Zhang (Slyz & Zhang, 2005) propose 2 compressors (BH and LM) for hyperspectral images. BH predicts a block as a scalar times the same block in the previous band. Coding contexts are defined by the quantized average error. LM predicts a pixel by choosing among different intraband predictions the one that works best for several pixels at the same position in previous bands.

Mielikainen and Toivanen proposed C-DPCM (Mielikainen & Toivanen, 2003), a method that classifies the pixels at the same location and through all the bands, with vector quantization. Interband prediction is performed using the pixels at the same position in 20 previous bands. Weights, calculated for each class/ band, are sent into the code, as well as the 2D template with the classes. Aiazzi et al. (Aiazzi et al., 1999) classify the prediction context of every

**Data**: Bands *I*(*z*)

**for** *every pixel I*(*z*) *<sup>x</sup>*,*<sup>y</sup>* **do** *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> LUTable(*<sup>I</sup>*

LUTable(*I*

0

**2.2 LAIS-LUT**

2

4

6

entropy

8

10

12

**end**

**Result**: Prediction for band *z*: *P*(*z*)

The Wavelet Transform as a Classification

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ) <sup>←</sup> *<sup>I</sup>*

Fig. 1. Look-up Table algorithm.

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* );

(*z*) *x*,*y*;

Criterion Applied to Improve Compression of Hyperspectral Images

Table 1. Average entropies for first scene of Jasper Ridge.

, *I*(*z*−1) and table initialized as LUTable(*i*) = *i*

**Original Wavelet S+P LUT prediction** 8.6656 6.6587 4.9504

> Original wavelet LUT

529

0 50 100 150 200 250

band

Fig. 2. Entropy values for every band of the Jasper Ridge image computed over the original image (dashed), over the prediction differences for the LUT method (gray), and over the 2D

An improvement over the LUT algorithm has been presented in (Huang & Sriraja, 2006). It was named LAIS-LUT after Locally Averaged Interband Scaling LUT and it behaves more accurately in presence of outliers. This modification adds an extra LUT table and the predictor

S+P wavelet transform (dotted). See averaged values in Table 1.

pixel using fuzzy clustering, and then calculate the weights for each class. For compression of hyperspectral images (Aiazzi et al., 2007) they divide each band into small blocks (16 × 16), they calculate weights for interband prediction over each pixel in a block, and then make a fuzzy clustering of the weights obtained for all the blocks. For each pixel a membership degree is computed according to the efficiency of the weights (from different clusters) on a causal context of the pixel. The final prediction for a pixel is obtained by a combination of linear predictions involving weights from different clusters, pondered by the degrees of membership of the pixel to each cluster. It is worth mentioning that wavelet-based compressors such as JPEG2000 (Taubman & Marcellin, 2002) have been successfully used for lossy compression of multiband images, either hyperspectral (Blanes & Serra-Sagrista, 2010; Fowler & Rucker, 2007) or general earth data (Kulkarni et al., 2006).

In this work we will improve the well-known algorithm which was developed for hyperspectral images: LAIS-LUT. This algorithm makes predictions of each pixel using other pixels in the same band. Which pixel is used for prediction is determined by inspecting in the previously encoded band. This algorithm uses two possible candidates for prediction. We will introduce the wavelet transform as a tool for classification in order to make better decisions about which of these two possible candidates acts as a more appropriate prediction. First, LAIS-LUT will be introduced, as well as the LUT algorithm on which it is based. Then we show how the wavelet transform is used to improve it. Finally, we give some results and conclusions of our method, called Enhanced LAIS-LUT.

#### **2. Look-up table based algorithms**

#### **2.1 LUT algorithm**

The Look Up Table algorithm (Mielikainen, 2006) is a fast compression technique based on predicting a pixel with another pixel from the same band. The previous band is inspected in order to determine which pixel in the same band is used for prediction. As an example, suppose pixel *I* (*z*) *<sup>x</sup>*,*<sup>y</sup>* in band *<sup>z</sup>* wants to be predicted. Then we seek on band *<sup>z</sup>* <sup>−</sup> 1 the pixel with the same intensity as *I* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* which is nearest to it in a causal neighborhood. Let *<sup>I</sup>* (*z*−1) *x*� ,*y*� be that pixel. Then, the prediction for pixel *I* (*z*) *<sup>x</sup>*,*<sup>y</sup>* will be *<sup>I</sup>* (*z*) *x*� ,*y*� . If no match is found, pixel *I* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* is the one selected for prediction. In order to speed things up, a look up table data structure is used for searching the pixel on the previous band. With this data structure the algorithm is efficiently implemented as shown in Fig 1 for consecutive bands *z* and *z* − 1. Then, the difference *<sup>I</sup>*(*z*) <sup>−</sup> *<sup>P</sup>*(*z*) between the band and its prediction is entropy coded and this process is repeated for *z* = 2, . . . , 224.

In Fig. 2 entropy values are ploted for each band of the Jasper Ridge image. In dashed line, entropies of the pixels of the image are ploted. 6 steps of the 2D S+P wavelet transform were computed and the entropy of the coefficients is plotted in dotted line for each band. Finally it is shown in gray line the entropy of the prediction differences for the LUT algorithm. It is remarkable how high is the compression achieved with this simple algorithm, which is only based on indexing and updating a table. It is entirely based on the premise of high correlation between bands and designed in order to take advantage of this fact.

**Data**: Bands *I*(*z*) , *I*(*z*−1) and table initialized as LUTable(*i*) = *i* **Result**: Prediction for band *z*: *P*(*z*) **for** *every pixel I*(*z*) *<sup>x</sup>*,*<sup>y</sup>* **do** *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> LUTable(*<sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ); LUTable(*I* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ) <sup>←</sup> *<sup>I</sup>* (*z*) *x*,*y*; **end**

Fig. 1. Look-up Table algorithm.

2 Will-be-set-by-IN-TECH

pixel using fuzzy clustering, and then calculate the weights for each class. For compression of hyperspectral images (Aiazzi et al., 2007) they divide each band into small blocks (16 × 16), they calculate weights for interband prediction over each pixel in a block, and then make a fuzzy clustering of the weights obtained for all the blocks. For each pixel a membership degree is computed according to the efficiency of the weights (from different clusters) on a causal context of the pixel. The final prediction for a pixel is obtained by a combination of linear predictions involving weights from different clusters, pondered by the degrees of membership of the pixel to each cluster. It is worth mentioning that wavelet-based compressors such as JPEG2000 (Taubman & Marcellin, 2002) have been successfully used for lossy compression of multiband images, either hyperspectral (Blanes & Serra-Sagrista, 2010; Fowler & Rucker,

In this work we will improve the well-known algorithm which was developed for hyperspectral images: LAIS-LUT. This algorithm makes predictions of each pixel using other pixels in the same band. Which pixel is used for prediction is determined by inspecting in the previously encoded band. This algorithm uses two possible candidates for prediction. We will introduce the wavelet transform as a tool for classification in order to make better decisions about which of these two possible candidates acts as a more appropriate prediction. First, LAIS-LUT will be introduced, as well as the LUT algorithm on which it is based. Then we show how the wavelet transform is used to improve it. Finally, we give some results and

The Look Up Table algorithm (Mielikainen, 2006) is a fast compression technique based on predicting a pixel with another pixel from the same band. The previous band is inspected in order to determine which pixel in the same band is used for prediction. As an example,

(*z*) *<sup>x</sup>*,*<sup>y</sup>* will be *<sup>I</sup>*

is the one selected for prediction. In order to speed things up, a look up table data structure is used for searching the pixel on the previous band. With this data structure the algorithm is efficiently implemented as shown in Fig 1 for consecutive bands *z* and *z* − 1. Then, the difference *<sup>I</sup>*(*z*) <sup>−</sup> *<sup>P</sup>*(*z*) between the band and its prediction is entropy coded and this process is

In Fig. 2 entropy values are ploted for each band of the Jasper Ridge image. In dashed line, entropies of the pixels of the image are ploted. 6 steps of the 2D S+P wavelet transform were computed and the entropy of the coefficients is plotted in dotted line for each band. Finally it is shown in gray line the entropy of the prediction differences for the LUT algorithm. It is remarkable how high is the compression achieved with this simple algorithm, which is only based on indexing and updating a table. It is entirely based on the premise of high correlation

between bands and designed in order to take advantage of this fact.

(*z*) *<sup>x</sup>*,*<sup>y</sup>* in band *<sup>z</sup>* wants to be predicted. Then we seek on band *<sup>z</sup>* <sup>−</sup> 1 the pixel

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* which is nearest to it in a causal neighborhood. Let *<sup>I</sup>*

(*z*) *x*�

,*y*� . If no match is found, pixel *I*

(*z*−1) *x*� ,*y*� be

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>*

2007) or general earth data (Kulkarni et al., 2006).

conclusions of our method, called Enhanced LAIS-LUT.

**2. Look-up table based algorithms**

that pixel. Then, the prediction for pixel *I*

**2.1 LUT algorithm**

suppose pixel *I*

with the same intensity as *I*

repeated for *z* = 2, . . . , 224.

**Original Wavelet S+P LUT prediction** 8.6656 6.6587 4.9504

Table 1. Average entropies for first scene of Jasper Ridge.

Fig. 2. Entropy values for every band of the Jasper Ridge image computed over the original image (dashed), over the prediction differences for the LUT method (gray), and over the 2D S+P wavelet transform (dotted). See averaged values in Table 1.

#### **2.2 LAIS-LUT**

An improvement over the LUT algorithm has been presented in (Huang & Sriraja, 2006). It was named LAIS-LUT after Locally Averaged Interband Scaling LUT and it behaves more accurately in presence of outliers. This modification adds an extra LUT table and the predictor is selected from one of the two LUTs. Using a scaling factor *α* which is precomputed on a causal neighbourhood, an estimate *<sup>P</sup>*˜(*z*) *<sup>x</sup>*,*<sup>y</sup>* is calculated for a current pixel *<sup>I</sup>* (*z*) *<sup>x</sup>*,*<sup>y</sup>* as *<sup>P</sup>*˜(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>=</sup> *<sup>α</sup><sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* , where

$$\alpha = \frac{1}{3} \left( \frac{I\_{i-1,j}^{(z)}}{I\_{i-1,j}^{(z-1)}} + \frac{I\_{i,j-1}^{(z)}}{I\_{i,j-1}^{(z-1)}} + \frac{I\_{i-1,j-1}^{(z)}}{I\_{i-1,j-1}^{(z-1)}} \right) \tag{1}$$

**3. The wavelet transform as a tool for classification**

Criterion Applied to Improve Compression of Hyperspectral Images

The Wavelet Transform as a Classification

work: orthogonal wavelets and lifting-based wavelets.

<sup>2</sup> <sup>−</sup><sup>1</sup> and detail coefficients **<sup>d</sup>** = [*sn*]

4 √2

prediction.

[*sn*]

*n*=0,..., *<sup>N</sup>*

filter is [*h*3, *h*2, *h*1, *h*0] = <sup>1</sup>

16), was also used in this work.

work are introduced.

**3.1 The wavelet transform**

In order to improve compression bitrates, the scaling factor in LAIS-LUT will be estimated more accurately. For that, a classification stage will be added to the algorithm. The scaling factor is computed as an average of quotients of collocated pixels in consecutive bands according to Equation 1. These pixels belong to a close neighbor of the pixel to be predicted (either on the current band, or on the previous band). Once classes are established for each pixel, not all the pixels in the close causal neighbor will be used for estimating the scaling factor *α*. Instead, we may use only those pixels that belong to the same class of the pixel to be predicted so as to obtain a more accurate estimation. This will enhance the LAIS-LUT algorithm, by allowing a better decision on which of the two look up tables yields a better

In order to establish classes we will make use of the wavelet transform. Since hyperspectral images have a considerable number of bands, the wavelet transform can be applied in the spectral direction. AVIRIS images have 224 bands. Considering each pixel as a vector in **Z**224, each of them may be 1D-wavelet transformed. And with the information of the wavelet transform of each 'pixel' classes can be determined. First, the wavelet transforms used in this

The wavelet transform allows the representation of a signal in multiresolution spaces. In the wavelet representation, the transformed signal can be viewed as an approximation of the original signal plus a sum of details at different levels of resolution. Each of these details and approximations are associated to function basis which have good time-frequency localization (Mallat, 1999). In images –a simple extension of the 1-dimensional case–, decorrelation is achieved obtaining a sparse representation. Two different types have been considered in this

For the classical orthogonal wavelet, consider a 1D signal **<sup>x</sup>** = [*xn*]*n*=0,...,*N*−<sup>1</sup> of length N (even). For a step of the wavelet, **x** is transformed into approximation coefficients **s** =

**x** with a lowpass filter followed by a decimation. The same process happens to **d**, but a high pass filter is used instead. For more steps, the wavelet transform is applied recursively over the approximation coefficients **s**. The original signal can be recovered if the inverse process is carried out in order (upsampling followed by the convolution with the corresponding reversed filter –see Fig. 4). Depending on the filter, wavelets with different properties can be obtained. In this work we use the Daubechies 4 wavelet (Daubechies, 1992) whose lowpass

*g*� = [−*h*0, *h*1, −*h*2, *h*3]. The Symmlet wavelet, with 8 vanishing moments (and filter of length

Another way of constructing different wavelets is by the lifting scheme (Daubechies & Sweldens, 1998). They are built in the spatial domain. The basic idea is to split the signal into two components, for instance, odd samples (**x**odd = *x*2*k*<sup>+</sup>1) and even samples (**x**even = *x*2*k*). Then, detail coefficients are obtained by using one component to predict the other: **d** = **x**odd − *P*(**x**even). Better predictions will yield more zeros, and therefore more decorrelation is

*n*=0,..., *<sup>N</sup>*

<sup>2</sup> <sup>−</sup>1, where **<sup>s</sup>** is the result of convolving

531

[<sup>3</sup> <sup>+</sup> *<sup>e</sup>*, 1 <sup>−</sup> *<sup>e</sup>*, 3 <sup>−</sup> *<sup>e</sup>*, 1 <sup>+</sup> *<sup>e</sup>*] with *<sup>e</sup>* <sup>=</sup> <sup>√</sup>3, and high pass filter

Since two values are now possible candidates for prediction (one for each LUT), the one that is closer to *<sup>P</sup>*˜(*z*) *<sup>x</sup>*,*<sup>y</sup>* is selected as the final prediction (Fig. 3 shows LAIS-LUT algorithm). When prediction *<sup>P</sup>*(*z*) for band *<sup>z</sup>* is estimated, the prediction error *<sup>I</sup>*(*z*) <sup>−</sup> *<sup>P</sup>*(*z*) is entropy coded. This is repeated for *<sup>z</sup>* <sup>=</sup> 2, . . . , 224. Notice that when the tables are not initialized, *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* is selected for prediction, and, the value of the first LUT is used for prediction only when this table is initialized.

**Data**: Bands *I*(*z*) , *I*(*z*−1) , tables LUTable1 and LUTable2 are not initialized **Result**: Prediction for band *z*: *P*(*z*)

**for** *every pixel I*(*z*) *<sup>x</sup>*,*<sup>y</sup>* **do**

**if** *both tables are not initialized on entry (x,y)* **then** *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> *<sup>α</sup><sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ; **else if** *only LUTable*<sup>1</sup> *is initialized in (x,y)* **then** *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> LUTable1(*<sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ); **else** *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> closer value of {LUTable1(*<sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ),LUTable2(*<sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* )} to *<sup>α</sup><sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ; **end** LUTable2(*I* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ) <sup>←</sup> LUTable1(*<sup>I</sup>* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ); LUTable1(*I* (*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ) <sup>←</sup> *<sup>I</sup>* (*z*) *x*,*y* ; **end**

Fig. 3. LAIS-LUT algorithm.

#### **2.3 LAIS-QLUT**

Mielikainen and Toivanen proposed a modification to the LAIS-LUT method in order to shrink Lookup tables (Mielikainen & Toivanen, 2008). For that, the value used for indexing into the LUTs is quantized and smaller LUTs can be used (i.e., the value *x* used as an index in the LUT is replaced by �*x*/*q*� for some quantization step *q*, being �·� : **R** → **Z** a function that maps a real number to a close integer). In the previous section, LAIS-LUT found exact matches in the previous band in order to obtain the predictor in the current band. For LAIS-QLUT, this search is no more exact and the predictor's selection in the current band is based on 'similarities' in the previous band. The best quantization step is obtained by an exhaustive search on each band (LAIS-QLUT-OPT version) or determined offline for each band by training on a set of images (LAIS-QLUT-HEU version). We decided not to add classes to the context over LAIS-QLUT, since a combination of the two would result in a prohibitive increase in complexity.

## **3. The wavelet transform as a tool for classification**

In order to improve compression bitrates, the scaling factor in LAIS-LUT will be estimated more accurately. For that, a classification stage will be added to the algorithm. The scaling factor is computed as an average of quotients of collocated pixels in consecutive bands according to Equation 1. These pixels belong to a close neighbor of the pixel to be predicted (either on the current band, or on the previous band). Once classes are established for each pixel, not all the pixels in the close causal neighbor will be used for estimating the scaling factor *α*. Instead, we may use only those pixels that belong to the same class of the pixel to be predicted so as to obtain a more accurate estimation. This will enhance the LAIS-LUT algorithm, by allowing a better decision on which of the two look up tables yields a better prediction.

In order to establish classes we will make use of the wavelet transform. Since hyperspectral images have a considerable number of bands, the wavelet transform can be applied in the spectral direction. AVIRIS images have 224 bands. Considering each pixel as a vector in **Z**224, each of them may be 1D-wavelet transformed. And with the information of the wavelet transform of each 'pixel' classes can be determined. First, the wavelet transforms used in this work are introduced.

#### **3.1 The wavelet transform**

4 Will-be-set-by-IN-TECH

is selected from one of the two LUTs. Using a scaling factor *α* which is precomputed on a

+ *I* (*z*) *i*,*j*−1 *I* (*z*−1) *i*,*j*−1

Since two values are now possible candidates for prediction (one for each LUT), the one that is closer to *<sup>P</sup>*˜(*z*) *<sup>x</sup>*,*<sup>y</sup>* is selected as the final prediction (Fig. 3 shows LAIS-LUT algorithm). When prediction *<sup>P</sup>*(*z*) for band *<sup>z</sup>* is estimated, the prediction error *<sup>I</sup>*(*z*) <sup>−</sup> *<sup>P</sup>*(*z*) is entropy coded. This is repeated for *<sup>z</sup>* <sup>=</sup> 2, . . . , 224. Notice that when the tables are not initialized, *<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* is selected for prediction, and, the value of the first LUT is used for prediction only when this table is

, tables LUTable1 and LUTable2 are not initialized

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ),LUTable2(*<sup>I</sup>*

Mielikainen and Toivanen proposed a modification to the LAIS-LUT method in order to shrink Lookup tables (Mielikainen & Toivanen, 2008). For that, the value used for indexing into the LUTs is quantized and smaller LUTs can be used (i.e., the value *x* used as an index in the LUT is replaced by �*x*/*q*� for some quantization step *q*, being �·� : **R** → **Z** a function that maps a real number to a close integer). In the previous section, LAIS-LUT found exact matches in the previous band in order to obtain the predictor in the current band. For LAIS-QLUT, this search is no more exact and the predictor's selection in the current band is based on 'similarities' in the previous band. The best quantization step is obtained by an exhaustive search on each band (LAIS-QLUT-OPT version) or determined offline for each band by training on a set of images (LAIS-QLUT-HEU version). We decided not to add classes to the context over LAIS-QLUT, since a combination of the two would result in a prohibitive increase in

+ *I* (*z*) *i*−1,*j*−1

*I* (*z*−1) *i*−1,*j*−1 ⎞

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* )} to *<sup>α</sup><sup>I</sup>*

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ;

(*z*) *<sup>x</sup>*,*<sup>y</sup>* as *<sup>P</sup>*˜(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>=</sup> *<sup>α</sup><sup>I</sup>*

⎠ (1)

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ,

causal neighbourhood, an estimate *<sup>P</sup>*˜(*z*) *<sup>x</sup>*,*<sup>y</sup>* is calculated for a current pixel *<sup>I</sup>*

⎛ ⎝ *I* (*z*) *i*−1,*j I* (*z*−1) *i*−1,*j*

*<sup>α</sup>* <sup>=</sup> <sup>1</sup> 3

where

initialized.

**else**

**end**

**end**

LUTable2(*I*

LUTable1(*I*

**2.3 LAIS-QLUT**

complexity.

Fig. 3. LAIS-LUT algorithm.

**Data**: Bands *I*(*z*)

**for** *every pixel I*(*z*) *<sup>x</sup>*,*<sup>y</sup>* **do**

*<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> *<sup>α</sup><sup>I</sup>*

, *I*(*z*−1)

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ;

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ) <sup>←</sup> *<sup>I</sup>*

**if** *both tables are not initialized on entry (x,y)* **then**

**else if** *only LUTable*<sup>1</sup> *is initialized in (x,y)* **then**

*<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> closer value of {LUTable1(*<sup>I</sup>*

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* ) <sup>←</sup> LUTable1(*<sup>I</sup>*

(*z*) *x*,*y* ;

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* );

(*z*−1) *<sup>x</sup>*,*<sup>y</sup>* );

**Result**: Prediction for band *z*: *P*(*z*)

*<sup>P</sup>*(*z*) *<sup>x</sup>*,*<sup>y</sup>* <sup>←</sup> LUTable1(*<sup>I</sup>*

The wavelet transform allows the representation of a signal in multiresolution spaces. In the wavelet representation, the transformed signal can be viewed as an approximation of the original signal plus a sum of details at different levels of resolution. Each of these details and approximations are associated to function basis which have good time-frequency localization (Mallat, 1999). In images –a simple extension of the 1-dimensional case–, decorrelation is achieved obtaining a sparse representation. Two different types have been considered in this work: orthogonal wavelets and lifting-based wavelets.

For the classical orthogonal wavelet, consider a 1D signal **<sup>x</sup>** = [*xn*]*n*=0,...,*N*−<sup>1</sup> of length N (even). For a step of the wavelet, **x** is transformed into approximation coefficients **s** = [*sn*] *n*=0,..., *<sup>N</sup>* <sup>2</sup> <sup>−</sup><sup>1</sup> and detail coefficients **<sup>d</sup>** = [*sn*] *n*=0,..., *<sup>N</sup>* <sup>2</sup> <sup>−</sup>1, where **<sup>s</sup>** is the result of convolving **x** with a lowpass filter followed by a decimation. The same process happens to **d**, but a high pass filter is used instead. For more steps, the wavelet transform is applied recursively over the approximation coefficients **s**. The original signal can be recovered if the inverse process is carried out in order (upsampling followed by the convolution with the corresponding reversed filter –see Fig. 4). Depending on the filter, wavelets with different properties can be obtained. In this work we use the Daubechies 4 wavelet (Daubechies, 1992) whose lowpass filter is [*h*3, *h*2, *h*1, *h*0] = <sup>1</sup> 4 √2 [<sup>3</sup> <sup>+</sup> *<sup>e</sup>*, 1 <sup>−</sup> *<sup>e</sup>*, 3 <sup>−</sup> *<sup>e</sup>*, 1 <sup>+</sup> *<sup>e</sup>*] with *<sup>e</sup>* <sup>=</sup> <sup>√</sup>3, and high pass filter *g*� = [−*h*0, *h*1, −*h*2, *h*3]. The Symmlet wavelet, with 8 vanishing moments (and filter of length 16), was also used in this work.

Another way of constructing different wavelets is by the lifting scheme (Daubechies & Sweldens, 1998). They are built in the spatial domain. The basic idea is to split the signal into two components, for instance, odd samples (**x**odd = *x*2*k*<sup>+</sup>1) and even samples (**x**even = *x*2*k*). Then, detail coefficients are obtained by using one component to predict the other: **d** = **x**odd − *P*(**x**even). Better predictions will yield more zeros, and therefore more decorrelation is

http://aviris.jpl.nasa.gov), we aim at determining classes for pixels. These images have 224 bands, each corresponding to a response in a certain range of the electromagnetic spectrum. Each pixel represents 20 meters and is allocated in a 2 byte signed integer. AVIRIS images usually have 614 columns and every 512 rows, the image is partitioned into 'scenes',

Considering each pixel *Ix*,*<sup>y</sup>* as a vector in **Z**<sup>224</sup> where each component belongs to a different band, different behaviours depending on what type of soil is being considered can be observed. Figure 6 shows the first scene of the Jasper Ridge image with four pixels belonging to different classes, marked with symbols '⊕', '◦', '♥' and '♦'. The spectral vectors associated to the same 4 positions are plotted in Figure 7; where the characteristic signal -called the

⊕

♥

533

°

In order to classify the image, for every pixel, a 1-D wavelet transform is applied along the spectral direction. The entropy of each transformed spectral vector is computed, giving an image *C* where *Cx*,*<sup>y</sup>* := *H*(*W*(*Ix*,*y*)), being *W*(·) a 1-D wavelet transform function and *H*(·) the entropy function weighted by wavelet subband. Since AVIRIS images have 224 bands, 5 steps of the wavelet transform are applied to each pixel. So we have 5 detail subbands and 1 approximation subband. Since decimation is performed after each step, these subbands are different in size. Therefore, the total entropy of the 1-D wavelet transform is computed as the

and each scene is stored in a different file.

The Wavelet Transform as a Classification

spectral signature- of the type of soil can be observed.

Criterion Applied to Improve Compression of Hyperspectral Images

♦

Fig. 6. Band 50 of Jasper Rigde image.

Fig. 4. Orthogonal wavelet analysis via convolutions and decimations followed by the synthesis via upsampling and convolution. Highpass filter *g* is obtained from lowpass filter *h* as *gk* = (−1)*kh*<sup>1</sup>−*k*, and *<sup>g</sup>*� indicates the *<sup>g</sup>* filter reversed.

achieved. In order to obtain the approximation part of the transformed signal, the unpredicted component is softened with an 'update' of the detail previously obtained: **s** = **x**even + *U*(**d**). See Fig. 5 for an illustration of this scheme. With the lifting scheme it is possible to construct wavelets that map integers into integers (Calderbank et al., 1998) by the use of a rounding operator at the prediction or update stage. For this work the (2, 2) and S+P wavelets have been considered (see Table 2).

Fig. 5. Lifting scheme analysis for a 1-D signal **x**.


Table 2. Wavelet transforms that maps integers into integers.

#### **4. Classification of hyperspectral images**

Using images captured by the scanner system called AVIRIS (Airborne Visible/Infrared Imaging Spectrometer) developed by JPL (Jet Propulsion Laboratory; see 6 Will-be-set-by-IN-TECH

**s**

**x**

been considered (see Table 2).

**x**

✲

*g*�

as *gk* = (−1)*kh*<sup>1</sup>−*k*, and *<sup>g</sup>*� indicates the *<sup>g</sup>* filter reversed.

✲

✓

✓

✒

✏<sup>↓</sup> <sup>2</sup> ✑

Fig. 4. Orthogonal wavelet analysis via convolutions and decimations followed by the synthesis via upsampling and convolution. Highpass filter *g* is obtained from lowpass filter *h*

**d**

achieved. In order to obtain the approximation part of the transformed signal, the unpredicted component is softened with an 'update' of the detail previously obtained: **s** = **x**even + *U*(**d**). See Fig. 5 for an illustration of this scheme. With the lifting scheme it is possible to construct wavelets that map integers into integers (Calderbank et al., 1998) by the use of a rounding operator at the prediction or update stage. For this work the (2, 2) and S+P wavelets have

✲

<sup>2</sup> (*x*2*<sup>n</sup>* <sup>+</sup> *<sup>x</sup>*2*n*+2) + <sup>1</sup>

<sup>4</sup> (*dn*−<sup>1</sup> <sup>+</sup> *dn*) + <sup>1</sup>

<sup>8</sup> (*sn*−<sup>1</sup> <sup>−</sup> *sn*) + <sup>3</sup>

Using images captured by the scanner system called AVIRIS (Airborne Visible/Infrared Imaging Spectrometer) developed by JPL (Jet Propulsion Laboratory; see

✒✑

2 

<sup>8</sup> (*sn* <sup>−</sup> *sn*+1) + <sup>2</sup>

8 *d* (1) *<sup>n</sup>*+<sup>1</sup> <sup>+</sup> <sup>1</sup> 2 

2 

✓✏

❄

❄

+ **s**

✲

✲

✲

✓

✓

✒

✏<sup>↑</sup> <sup>2</sup> ✑

*g*

✻

**d**

✒✑

✓✏<sup>+</sup> ✲ **<sup>x</sup>**

❄

*h*

✲

✏<sup>↑</sup> <sup>2</sup> ✑

✲

✒

✲

✲

✏<sup>↓</sup> <sup>2</sup> ✑

✲

✒

✲

**Wavelet Formula**

 1

 *d* (1) *n* 2 

 1

**x**[2*k* + 1]

**x**[2*k*]

Fig. 5. Lifting scheme analysis for a 1-D signal **x**.

(2, 2) *dn* = *x*2*n*+<sup>1</sup> −

S+P *d*

*sn* = *x*2*<sup>n</sup>* +

*sn* = *x*2*<sup>n</sup>* +

*dn* = *d*

**4. Classification of hyperspectral images**

(1) *<sup>n</sup>* <sup>=</sup> *<sup>x</sup>*2*n*+<sup>1</sup> <sup>−</sup> *<sup>x</sup>*2*<sup>n</sup>*

(1) *<sup>n</sup>* + 2

Table 2. Wavelet transforms that maps integers into integers.

split *P U*

✻

✒✑ ✓✏–

✻

✲

*h*�

✲

http://aviris.jpl.nasa.gov), we aim at determining classes for pixels. These images have 224 bands, each corresponding to a response in a certain range of the electromagnetic spectrum. Each pixel represents 20 meters and is allocated in a 2 byte signed integer. AVIRIS images usually have 614 columns and every 512 rows, the image is partitioned into 'scenes', and each scene is stored in a different file.

Considering each pixel *Ix*,*<sup>y</sup>* as a vector in **Z**<sup>224</sup> where each component belongs to a different band, different behaviours depending on what type of soil is being considered can be observed. Figure 6 shows the first scene of the Jasper Ridge image with four pixels belonging to different classes, marked with symbols '⊕', '◦', '♥' and '♦'. The spectral vectors associated to the same 4 positions are plotted in Figure 7; where the characteristic signal -called the spectral signature- of the type of soil can be observed.

Fig. 6. Band 50 of Jasper Rigde image.

In order to classify the image, for every pixel, a 1-D wavelet transform is applied along the spectral direction. The entropy of each transformed spectral vector is computed, giving an image *C* where *Cx*,*<sup>y</sup>* := *H*(*W*(*Ix*,*y*)), being *W*(·) a 1-D wavelet transform function and *H*(·) the entropy function weighted by wavelet subband. Since AVIRIS images have 224 bands, 5 steps of the wavelet transform are applied to each pixel. So we have 5 detail subbands and 1 approximation subband. Since decimation is performed after each step, these subbands are different in size. Therefore, the total entropy of the 1-D wavelet transform is computed as the

Enhanced LAIS-LUT LAIS- LUT

535

Image Daubechies 4 Symmlet 8 S+P (2,2) LUT

Jasper Ridge 3.43895 3.43855 3.43914 3.43833 3.42 3.23 Cuprite 3.71393 3.71394 3.71390 3.71390 3.58 3.44 Lunar Lake 3.58201 3.58198 3.58198 3.58201 3.53 3.4 Table 3. Compression ratios for Enhanced LAIS-LUT (with classes for estimating the scaling

sum of the entropies of each subband, weighted by the size of the subband relative to the size

In Figure 8 a grayscale image of the entropies is displayed. This image may be further split into classes with an unsupervised classifier. In this work, mean-shift (Comaniciu & Meer,

Results of compressing AVIRIS images with LAIS-LUT and the enhanced version (in which the scaling factor is calculated with pixels belonging to the same class, considering a causal

The name of the wavelet in the table indicates the wavelet used to transform each pixel (and all the bands) followed by the entropy estimation and mean shift classification. Compression ratio results for LAIS-LUT and LUT algorithms were obtained from (Mielikainen & Toivanen, 2008). It can be observed that the enhanced version of LAIS-LUT using Daubechies 4 for

We may conclude that the scaling factor *α* plays an important role in the compression performance of LAIS-LUT algorithm. When introduced, it was intended to decrease the deterioration produced by outliers in the original LUT algorithm. We have also been able to make use of the information in the wavelet domain and apply it to develop an efficient classifier. Since hyperspectral images have many bands because of their high spectral resolution, the information of the signal that each pixel represents (in all bands) was well captured by the wavelet transform and was fed into a powerful classifier such as mean-shift,

Aiazzi, B., Alba, P., Alparone, L. & Baronti, S. (1999). Lossless compression of

Aiazzi, B., Alparone, L., Baronti, S. & Lastri, C. (2007). Crisp and fuzzy adaptive spectral

Blanes, I. & Serra-Sagrista, J. (2010). Cost and scalability improvements to the karhunen-loêve

multi/hyperspectral imagery based on a 3-D fuzzy prediction, *IEEE Trans. Geos.*

predictions for lossless and near-lossless compression of hyperspectral imagery, *IEEE*

transform for remote-sensing image coding, *Geoscience and Remote Sensing, IEEE*

neighbor of 4 × 4 pixels around the one to be predicted) is shown in Table 3.

Low Altitude 3.17875 3.17546 3.17793 3.17797

factor) and classical LAIS-LUT and LUT algorithms.

Criterion Applied to Improve Compression of Hyperspectral Images

The Wavelet Transform as a Classification

classification outperforms the other methods.

*Remote Sensing* 37(5): 2287–2294.

*Transactions on* 48(7): 2854 –2863.

*Geos. Remote Sensing Letters* 4(4): 532–536.

giving good compression results.

**6. References**

of the whole transform.

2002) was used.

**5. Results**

Fig. 7. Spectral signatures of different classes: values of selected pixels in Fig. 6 throughout the bands.

Fig. 8. Entropy of the wavelet transform of every pixel. Darker pixels indicate lower entropy values.


Table 3. Compression ratios for Enhanced LAIS-LUT (with classes for estimating the scaling factor) and classical LAIS-LUT and LUT algorithms.

sum of the entropies of each subband, weighted by the size of the subband relative to the size of the whole transform.

In Figure 8 a grayscale image of the entropies is displayed. This image may be further split into classes with an unsupervised classifier. In this work, mean-shift (Comaniciu & Meer, 2002) was used.
