SciELO - Scientific Electronic Library Online

 
vol.43 número3Set-membership estimation theory for coupled mimo Wiener-like modelsHomography-based pose estimation to guide a miniature helicopter during 3d-trajectory tracking índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Latin American applied research

versión On-line ISSN 1851-8796

Lat. Am. appl. res. vol.43 no.3 Bahía Blanca jul. 2013

 

Correlation-based inter and intra-band predictions for lossless compression of multispectral images

D. Acevedo and A. Ruedin

Departamento de Computación, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires
C1428EGA Buenos Aires, Argentina. Email: dacevedo@dc.uba.ar, ana.ruedin@dc.uba.ar

Abstract — We present a new lossless compressor for multispectral images having few bands. The mentioned compressor takes into account variations in spectral correlation in order to determine the appropriate spectral and spatial prediction to be performed. The algorithm exploits 2 different facts. On one hand, highly correlated bands may be efficiently compressed with fast computations. On the other hand, a class-conditioned wavelet-based compressor, which is more time-consuming, has given very high compression ratios, even in the case of lowly correlated bands. Our correlation dependent hybrid algorithm yields high compression ratios, outperforming state-of-the-art lossless compressors, and has reasonable execution times.

Keywords — Prediction; Wavelet Coefficients; Classification; Multispectral; Lossless Compression.

I. INTRODUCTION

Satellites continually feeding images to their base, pose a challenge as to the design of compression techniques to store these huge data volumes. Even though the price of storage devices tends to drop, the cost of storing satellite images continues to be the dominant cost of the system. Proper use of compression can reduce this cost.

We aim at lossless compression of multispectral images having few bands, such as Landsat 7 ETM+ images. These images consist of several images (or bands) obtained by filtering radiation from the earth at different wavelengths. Compression is generally achieved thro-ugh reduction of spatial as well as spectral correlations.

To reduce correlations, prediction-based methods have been applied. For each pixel being encoded, a prediction for that pixel is performed, using information from a neighborhood of the pixel called the prediction context. The differences between the predictions and the original pixels are then encoded with an entropy based coder. The more accurate the predictions, the smaller the differences, and the higher the compression ratio. To exploit interband correlations, predictions involve pixels from other bands. In the final coding step, sometimes the coding context (another neighborhood) of the pixels is used: different coding contexts define different coders. Class information is also used: by conditioning the predictions to the class of the prediction context, or coding the prediction differences conditioned to the class of the coding context, compressors try to capture variations in the local statistics of the image. In order to perfectly recover the information at the decompressor's side, both the prediction context and the coding context must consist of pixels that have already been coded, and are available to the decoder. In the case of linear predictions, the weights for the prediction are calculated over a set of pixels called the training set. Then either the weights are written into the code, or they must be recalculated by the decoder: in the latter case the training set must also consist of already coded pixels.

JPEG-LS (Weinberger et al., 2000) uses a nonlinear median edge prediction, and a large number of coding contexts based on vertical and horizontal gradients. We propose an extension of the median edge predictor to multispectral images, and apply it to compress the highly correlated thermal bands (6.1 and 6.2) of Landsat 7 images (Ruedin and Acevedo, 2005).

Another approach to reduce correlations in a 2D image is through the application of wavelet transforms (Mallat, 1999). By representing an image as the sum of its details at different resolutions and orientations (plus a coarse approximation), a wavelet transform substantially reduces spatial correlations. Detail subbands have histograms that are peaked and centered at zero: their entropies are small. This is a remarkable property of the wavelet transform, and explains why wavelets are in the state of the art for image compression.

Traditional wavelets produce real coefficients. For lossless compression, wavelets that map integers into integers have been designed (Calderbank, 1998). These transforms are reversible when the values of the original image are integers, as is the case. They are used in SPIHT (Said and Pearlman, 1996), and JPEG 2000 (lossless) (Rabbani and Joshi, 2002).

Landsat images have 8 bands, and around 7000 rows and columns per band. In Table 1 we have interband correlations on a 256×256 sample, where we can observe that correlation values lie in a wide range(-0.07 to 0.98). Because Landsat images have few bands and low spectral resolution, and because the wavelengths at which the bands are filtered are irregular- see Table 2- the mentioned wavelet-based algorithms, when extended for volumetric data, such as 3D-SPIHT (Tang et al., 2003) or cube JPEG2000 (Lee et al.,2000), have given poor results on these images. DEC3 (Benazza, 2002), another algorithm, applies a 2D wavelet transform to each band: a linear prediction involving pixels from several previous bands, is built into the lifting scheme of the wavelet transform. Yet it does not use class information, which might improve results.

If correlation is high, pixels in the actual block are predicted with an extension of the median edge predictor. This extension includes pixels in the previous band, after equalization of one block to the other. If correlation with the previous block is low, spatial correlations in the block are reduced via an integer wavelet transform, and both spatial and spectral correlations are further reduced by prediction of wavelet coefficients. Finally prediction differences are encoded with a context-based arithmetic coder. By switching between the 2 algorithms we decrease the computational burden, whilst maintaining high compression ratios.

This paper is organized as follows. In Section II we explain the algorithm for block ordering. In Section III we address the extension of the median edge predictor. In Section IV we give the wavelet-based prediction algorithm, including the classification step, necessary for calculating weights tuned to the class. Numerical tests and conclusions are given in Section V and VI.

II. BLOCK ORDERING

Our compressor processes and encodes one stack after the other, in raster scan order. Inside a stack, to compress a given block, our compressor switches between 2 different algorithms, depending on the correlation between this block and the previous one. Both algorithms are based on predictions. More accurate predictions will give lower prediction differences and higher compression ratios. The order of the blocks affects the accuracy of predictions; it is therefore desirable to have an optimal ordering of the blocks, so that highly correlated blocks are coded consecutively. This ordering may not be unique and must be sent to the decoder (into the header of the stack) to recover all the blocks in the right order.

We apply a greedy heuristic to find an ordering of the blocks. Inside a stack of blocks, the mechanism is as follows: at every step the algorithm chooses the block having the greater correlation with the previous one. The algorithm is initialized with each block k in the stack, so that we have as many different orders as blocks there are in the stack. From the all the orders found, we select the one having the greatest sum of correlations.

III. CODING OF HIGHLY CORRELATED CO-LOCATED BLOCKS IN CONSECUTIVE BANDS

An algorithm was designed for the compression of blocks which were highly correlated to the previous block in the same stack. These blocks have exactly the same location but belong to different bands. It is based on histogram equalization and context-based median edge prediction (Weinberger et al., 2000).

Table 1 : Spectral correlations in a stack of blocks

Table 2: Wavelengths captured by each band.

Correlation measures the linear relationship of 2 random variables. In two highly correlated blocks, co-located pixels are linearly related; hence their histograms are similar. This is not true when correlations are low. Compare the histograms in Fig. 1 and those in Fig. 2.


Figure 1: Histograms of two co-located blocks in consecutive bands having correlation -0.06.


Figure 2: Histograms of two co-located blocks in consecutive bands having correlation 0.98.

The median edge predictor, applied to a single image, predicts a pixel by inspecting the 3 nearest pixels, and calculates the median between the pixel up, the pixel left, and value of the interpolating plane (passing through the 3 pixels), at position (i, j) (see Fig. 4):


Figure 4: Multispectral image: pixels used for predicting

Then the difference between the actual pixel and its prediction is encoded.

For the compression of a block highly correlated to the previous block, a predictor similar to the median edge predictor was sought, which would also include the information of pixel , that is, the pixel at the same position as the one to be predicted, in the previous block. This was insufficient. Histograms from block k and block k−1, however similar in shape, have different mean values and different standard deviations. It was necessary to equalize block k−1 to block k, by matching their cumulative histograms, before we made use of pixel .

Consider block k as a sample of a random variable U, with cumulative distribution FU, and consider block k−1 as a sample of another random variable V, with cumulative distribution FV. For each value u, let v be the value such that FV(v)= FU(u), and define . Then g(U) has the desired cumulative distribution, similar to FV (Gonzalez and Woods, 2002).

By calculating g(u) for each pixel u of block k−1, we transform block k−1 to a new image having a cumulative histogram similar to that of block k. In Fig. 3 we have the resulting histogram of equalizing the first block (whose original histogram is in Fig. 2(a) to the second block (whose original histogram is in Fig. 2(b)). The height of the histogram values remains the same as in Fig. 2(a), but the range of pixel values has changed to those of Fig. 2(b).


Figure 3: Histogram of first block equalized to second block (with reference to Fig. 2)

Now the proposed prediction of our method for pixel is the median of 4 values:

(1)

After applying formula (1), the predictions errors - were encoded with a context-based arithmetic coder. The contexts used were two horizontal gradients above the pixel to be predicted: and , plus one vertical gradient to the left: These contexts are 3 dimensional vectors of integer components. To decrease the great number of possible context classes, the values of the gradients were quantized into three groups , giving a total number of 27 context classes.

This method has low complexity; a block is coded (or decoded) in less than 9 milliseconds, and a whole image is compressed (or decompressed) in 44.42 seconds.

IV. CODING OF LOWLY CORRELATED CO-LOCATED BLOCKS IN CONSECUTIVE BANDS

For blocks having low correlations, we applied a different algorithm based on an integer wavelet transform and linear predictions, called CLWP (Class-conditioned Lossless Wavelet-based Predictive compressor) (Ruedin and Acevedo, 2010).

A. Integer Wavelet Transform on Blocks

We applied a wavelet transform that maps integers into integers to each block. For our compressor we chose to apply the S+P transform (Said and Pearlman, 1996), used in SPIHT.

In Fig. 5 we have (a) an original block, and (b) the S+P transformed block, in which every subband hasbeen rescaled. In Fig. 6(a) are the names of the wavelet subbands: e.g. V1, H1, D1 are respectively vertical, horizontal and diagonal fine detail coefficients. An increase in the indexes indicates less resolution. A3 consists of coarse approximation coefficients.

In our compressor, 5 steps of the S+P transform are computed on each block, giving 16 subbands. The order in which the subbands are encoded is the following: A5 (approximation coefficients), {Vk, Hk, Dk} for k =5, ..., 1 (detail coefficients). Each subband is encoded in raster scan order.


Figure 5: A 256×256 block and its S+P Wavelet transform (3 steps).


Figure 6: Wavelet transform, 3 steps (a) Wavelet sub-bands. (b) S+P Wavelet transform of original block, with classified fine detail coefficients (black=water).

Instead of encoding the wavelet coefficients directly, each wavelet coefficient is predicted with an affine combination of already coded coefficients. Prediction differences are entropy coded. The prediction is carried out only for the detail coefficients of the two finest scales, that represent (15/16)100= 93.75% of the coefficients. The remaining wavelet coefficients are encoded as they are.

Distributions of wavelet coefficients vary with their resolution level, their orientation, and the class, so that we perform different predictions for each band/ block/ subband/ class.

B. Classification

We have chosen a K-means clustering algorithm, which does not require previous knowledge on landscapes. For reasons of speed, we chose a classification of pixels (involving one block) rather than vectors (involving all the blocks in a stack), so that both classification and weights (for each landscape and subband) were calculated anew for each block. There was no significant increase in the compression rates when more than 2 classes were considered: accordingly, the pixels in each block were separated into 2 classes.

For reconstruction it is necessary to have exactly the same classes at the decoder: this required either sending the template of the pixel classes to the decoder (which would constitute a prohibitive overhead), or have both the coder and the decoder determine the class in an identical way. This could only be done if the classification involved data that was available to the decoder: therefore we based our prediction on the classification of the corresponding block in the band previous to the one being encoded.

That is, assuming correlation between 2 consecutive bands to be significant, we conjecture that the pixels at the same position in 2 consecutive bands, generally belong to the same class. This conjecture is also applied to initialize K-means in a block of a given band with the classes obtained in the corresponding block in the previous band. The decompressor is faster than the compressor, because it is given the weights for the prediction.

Once we have separated the pixels of the image block into 2 classes, we have a template indicating the class of each pixel. This induces a classification of the wavelet coefficients in the 6 predicted subbands. We first construct a second template that is a fourth (in size) of the original one, to indicate the classes in subbands V1, H1 and D1. For this, we divide the original template in 2×2 submatrices, and the most numerous class in a submatrix is an entry in the second template (see Fig. 7). Similarly, we obtain a 3rd template indicating the classes in subbands V2, H2 and D2.


Figure 7: Left: Original template. Right: Second template.

K-means was applied to separate the pixels of Fig. 5 (a) into land and water. With our procedure we obtained a classification of the fine detail coefficients in Fig. 5 (b). The result may be observed in Fig. 6 (b).

C. Prediction of Wavelet Coefficients

We have analyzed different models for the statistical dependencies between wavelet coefficients of a multispectal image. Our best results were found by modeling wavelet detail coefficients in a given subband and class, as an affine combination of (i) neighboring coefficients up and left (in the same class) and (ii) the wavelet coefficient from the previous spectral band, at the same location.

Weights for the affine prediction are calculated by least squares over all the coefficients in the same sub-band belonging to the same class. They are sent into the bitstream. We calculate the affine prediction for each wavelet detail coefficient, round the result and encode the prediction differences with a context-based arithmetic coder.

Each stack of blocks is coded independently. The encoder deals with 2 different types of values: i) the prediction differences for the wavelet coefficients belonging to the fine detail scales ii) the wavelet coefficients themselves, belonging to the coarser detail sub-bands and the approximation subband, for which prediction was not worthwhile. All values are coded with an adaptive arithmetic coder. For each block there are 16 wavelet subbands. Since the entire subbands are coded, the arithmetic coder has to be restarted 16 times per block. Prediction differences of the 6 fine detail subbands are coded conditioned to the class. CLWP needs 73 milliseconds to compress a block, and 6 minutes to compress a whole image (3 min.12 sec.to decompress it) on a PC running at 2.4 Ghz.

V. NUMERICAL RESULTS

A block diagram of our proposed method is given in Fig. 9. The decision threshold that decides which algorithm is executed for each block, based on inter-band correlation, is calculated empirically by training on several images.


Figure 9: Block diagram. AC stands for Arithmetic Coding.

Table 3 shows the results of testing our compressor on 13 Landsat 7 ETM+ images1, where P/R-AD stands for Path/Row-Adquisition Date (yyyy-mmdd). We compare the compression ratio of our proposed method to different lossless compressors, such as SLSQ-OPT (Rizzo et al., 2005), Differential JPEGLS (in which the difference between 2 consecutive bands are coded with JPEG-LS), JPEG 2000 Color, 3D SPIHT, DEC3, KLT+DWT (reversible Karhunen-Loêve transform is applied in the spectral dimension; then each band is 2D-transformed with an integer wavelet), LUT (Mielikainen, 2006) and LAIS-LUT (Huang and Sriraja, 2006) (these last two, developed for multiband hyperspectral images). Our proposed method outperforms the others. LUT and LAIS-LUT, which give the best results when applied to hyperspectral AVIRIS images, rely exclusively on interband prediction: this is why here they are not so effective, interband correlation being sometimes very low.

In Fig. 8 we have a Landsat image, where each square represents a stack of blocks. The shade of each square reveal show many blocks in the stack were processed with the first algorithm (extended median edge prediction), i.e. how many blocks in the stack were highly correlated to the co-located block in the previous band in the new ordering. This is for values inside the image; zero-valued pixels in the borders of the image were detected and encoded with run-length encoding. 43% of the blocks were coded with the extended median edge predictor. This represented a reduction in 37.5% of the time employed to code all the image with the wavelet- based predictor CLWP.


Figure 8: Each square represents a stack of blocks. The shade indicates the amount of blocks for which interband correlation was above the threshold.

1. We thank CONAE (Comisión Nacional de Actividades Espaciales) for providing the images.

Table 3: Compression ratios obtained with different lossless compressors.

VI. CONCLUSIONS

In this paper we have presented a new lossless compressor that exploits both spatial and spectral correlations present in multispectral images. Highly correlated blocks are coded with a simple pixel median prediction involving pixels from the previous band after equalization. Lowly correlated blocks are coded with algorithm CLWP: with an integer wavelet transform we obtain a significant reduction of spatial correlation. To reduce spectral correlation (and remaining spatial correlation) we rely on inter and intraband predictions and encoding of prediction differences. These are performed on the wavelet coefficients. Conditioning the estimation of weights to the class, band, block and sub-band, has led to more accurate predictions. This has resulted in lower entropies of prediction differences.

REFERENCES
1. Benazza-Benyahia A., J. Pesquet and M. Hamdi, "Vector-lifting schemes for lossless coding and progressive archival of multispectral images," IEEE Tr. Geo. Rem. Sens., 40, 2011-2024 (2002).
2. Calderbank, A.R., I. Daubechies, W. Sweldens and B. Yeo, "Wavelet transforms that map integer to integers," Applied and Computational Harmonics Analysis, 5, 332-369 (1998).
3. Gonzalez, R.C. and R.E. Woods, Digital Image Processing, Prentice Hall (2002).
4. Huang, B. and Y. Sriraja, "Lossless compression of hyperspectral imagery via lookup tables with predictor selection," Proceedings of SPIE, 6365 (2006).
5. Lee, H.S., N.-H.Younan and R.L. King, "Hyperspectral image cube compression combining JPEG 2000 and spectral decorrelation," Proc. IEEE IGARSS, 6, 3317-3319 (2000).
6. Mallat, S., A Wavelet Tour of Signal Processing, Academic Press (1999).
7. Mielikainen, J., "Lossless compression of hyperspectral images using lookup tables," IEEE Signal Processing Letters, 13, 157-160 (2006).
8. Rabbani, M. and R. Joshi, "An overview of the JPEG 2000 still image compression standard," Sig. Proc. Im. Comm., 17, 3-48 (2002).
9. Rizzo, F., B. Carpentieri, G. Motta and J.A. Storer, "Low complexity lossless compression of hyperspectral imagery via linear prediction," IEEE Signal Processing Letters, 12, 138-141 (2005).
10. Ruedin, A. and D. Acevedo, "Prediction of coefficients for lossless compression of multispectral images," Proceedings of SPIE, 5889, 202-211 (2005).
11. Ruedin, A. and D. Acevedo, "A class-conditioned lossless wavelet-based predictive multispectral image compressor," IEEE Geoscience and Remote Sensing Letters, 7, 166-170 (2010).
12. Said, A. and W. Pearlman, "A new fast and efficient image codec based on set partitioning in hierarchical trees," IEEE Trans. Circuits Sys. Video Tech., 6, 243-250 (1996).
13. Tang, X., S. Cho and W.A. Pearlman, "Comparison of 3D set partitioning methods in hyperspectral image compression featuring an improved 3D-SPIHT," Proc. Data Compression Conf. (2003).
14. Weinberger, M., G. Seroussi and G. Sapiro, "The LOCO-I lossless image compression algorithm: principles and standardization into JPEG-LS," IEEE Trans. on Image Processing, 9, 1309-1324 (2000).

Received: April 25, 2012
Accepted: October 10, 2012
Recommended by Subject Editor: Gastón Schlotthauer, María Eugenia Torres and José Luis Figueroa.