Study on Lossless Compression of Multispectral Remote Sensing Images
Kai Xu,Hiroshi Okayama
Remote Sensing and Image Research Center
Chiba University I-33 Yayoi-Cho, Chiba 260 Japan
Koji Kajiwara
Institute of Industrial Science, University of Tokyo 7-22 Roppongi
Minato-ku, Tokyo 106, Japan
Abstract
This paper describes a partially recoverable lossless data compression scheme for multispectral remote sensing images of raster format. The method takes both the intrepid and interurban redundancy into account and deals with intraband redundancy by converting each band of image to its differential image. The integrand correlations are handled by pre-coding the combination image (pairs of differential values) of two differential images. Then a line-by-line two bands of original images. The coding system is designed to maintain a sufficient line index information which ensures the possibility of reconstruction any part of the original image with no need to decode the whole compressed data ached. Thus we can access the compressed image data just as we do it with the original one. the method was applied to SPOT HRV and LANDSAT TM multispectral images an overall compression ratio up to 3.3. and 2.3 respectively was achieved.
Introduction
Multispectral satellite images, such as SPOT, LANDSAT TM or MSS images, are rich resources for earth observations and environmental analysis. However, satellite image data are normally of enormous sizes. For example, a single SPOT HRV multispectral (3 bands) scene takes up to a total 30 Mbytes of storage facilities. The steady growth implies a significant increase of storage requirements for such sensors, in the near future. from the point of view of image interpreters or image analysts who work on earth observation data, a global earth observation image data-base is highly desired. Thus, an easily accessible data format of digital images with low storage requirements should contribute somewhat to an efficient data-base management scheme.
For single band images, although a rather high sampling rate and quantization is required to obtain subjectively high-quality pictures, statistical analysis of images indicates that the entropy of most pictures is substantially less than what these high sampling and quantization rates requie (Tzay Y. Young, 1986). This implies that a rather high redundancy exists among the neighboring pixels due to the scanning mechanism of most sensors. For the multispectral remotely sensed images, interband redundancy should also exist as well because of the similar reflectance received from the objects for neighboring bands. And also, from an intuitive point of view, the fact that the objects in all bands share the same boundaries should lead to a considerable integrand correlations. These make the compression of multispectral images desirable and possible for storage reduction.
Generally, image compression schemes can be divided in two categories. One is the lossy approach by which a high compression ratio can be achieved when allowing a certain limited distortion in the reconstructed image provide that it leaves no apparent difficulty four visual interpretations. For this case the transform coding or a hybrid of DPCM (Differential Pulse Code Modulation) with transform coding scheme has been proved as practical and is widely utilized (Tzay Y. Young, 1986). The alternation is usually considered when a large amount of data is to be stored for future retrieval and analysis. For the lossless image compressions a reversible DPCM or its adaptive extension followed by a variable length encoder is widely used.
Our consideration for a lossless image compression concentrates on two objectives. One is to store the remotely sensed images in a lower storage volume for future retrieval. The other is to make the compressed file as easily accessible as the uncompressed original ones. Specifically, the compressed data must be partially recoverable so that we can extract any line, any block of nay band of original images with no need to extend the whole compressed files ahead. Thus, when the decoder is incorporated into some software developed in our center the users can just treat the compressed images the same as uncompressed ones. This is accomplished by maintaining a sufficient extra information as indices. Experimental results show that the overall ( including all the index tables ) compression ratio for SPOT HRV ( full scene 5369x3000 pixels/band x 3bands, 8 bits/pixel) and LANDSAT TM ( screen size 512x400.
Pixels/band x 4bands, 8bits/pixel) up to 3.3. respectively is achievable.
System design
1 Differential image
The most commonly used for of predictive coding system is DPCM, which uses a linear predictor to predict the incoming signal based on a weighted sum of the adjacent elements. The DPCM generates a set of uncorrelated signals, referred to as differential signals, from the correlated data, and the processing could either be 2-dimensional or I- dimensional. Since our goal is to make the compressed file partially decodable, we here choose the I-dimensional DPCM to reduce the redundancy degree among the neighboring pixels on each scan line. Thus given the image size (N pixels?line x M lines ) the inverse processor can easily reproduce any original data segment on any line. The prediction of pixel at ith column and jth line can be expressed like (1),
The optimal performance of the prediction ( in Mean Square Error sense ) depends on the optimal choose of the prediction coefficients, the aks in (tzay Y. Yung, 1986 ). Some statistical studies showed (Robert E Rice, 1971 ) that only the ui-lj contributed significantly to the prediction of ^u
ij. Thus (1) can be simplified to as
which is simply a differencing operation. We view the new image consists of the u
0js and the prediction error e
ijs as the differential image. Our next step is to code this differential image instead of the original, since the differential image, generally, has more compact dynamic range than that of the original image. Here a truncating operation is necessary to be carried out to the e
ijs before the coding procedure because, in most cases, remote sensing image data are usually quantized to an 8 bits/pixel level so that the value of e
ijs may vary from - 255 to + 255, which requires 9 bits for every e
ij, and this is not convenient for implementation on computers. We truncate the e
ijs according to (3) where the T(e
ij) is the truncated value of e
ij, and the table h
ij holds the transformed value of e
ijs which are out of the range [-126,+126] and the end-of-line (EOL) marks with value 129. The T (e
ij)s are to be recorded in a file as the actual differential image for further coding. And the h
ijs must be stored as a reference table for decoding.
So the inverse operation will simply be the inverse of T(e
ij) with a table looking up into hij whenever the value of 253 or 254 is encountered in the coded file. A line - by-line inverse of the truncating operation can be accomplished by counting the EOLs in h
ij.
2 Combination image
By converting the individual band image to its differential image the intraband correlations are significantly reduced. As explained above there should exist some correlations between adjacent bands, termed as the interband correlations. It is quite natural to consider a similar processing of interband prediction for corresponding pixels in two
Differential images then coding the prediction error image as proposed by (Soon Ing Yann, 1991) or simply the band-difference image. However, results of simulations showed that this did not work satisfactorily; Experiments indicated that simply taking the band-difference did not work well for correlation reduction. There must exist an optimal prediction coefficient for two band images, but to calculate it needs some intensive statistical computations especially when the data are of large volume
To reduce the correlations between two images we here try a method of pre-coding the combination value of corresponding pixel intensity values in two differential images. If the data file containing the combination values is viewed as an image file we call it the combination image. To guarantee the reversibility the 2-dimensional histogram of the underlying differential images must be stored. We allocate a 16 bits memory for each combination value in which the upper byte holds the pixel value of band one differential image and the lower byte for the corresponding pixel in
Band two. This can be written as (4) , where c
ij is the combination value for pixels at ith line jth column and uk
ij stands for the pixel intensity value of band k differential image.
Clearly, when the u
kijs have 256 possible values (gray levels ) the c
ijs will have 65536 possible values and to code the c
ijs directly will actually achieve nothing. Experiments showed that when (4) was applied to multispectral remote sensing images after the pixel differencing operation described above the number of different c
ij values was much less than 65536. This lead us to think of assigning a particular number to each different combination value was much less representation. A simple way of doing this is to first sort the 2-dimensional histogram calculated from the underlying differential images according to the decreasing of occurrence of each pixel combination and then use the index number of sort number to represent different combinations. For example, let L be the number of different pixel combinations so each of them can be represented by a integer in (0,1,2,…,L-1) where 0 stands for the most frequently occurred combination and 1 for the second and so on. By far a converting table is constructed and now we can scene the two differential images simultaneously and convert each pixel combination according to the converting table to generate a new data file. The reverse processing is also guaranteed as long as the sorted 2-D histogram of length L is stored.
Consider each element as a point in a plane with the X-Y coordinates be its lower byte and upper byte respectively then this new data file will be the points covering a rectangular area in the X-Y plane with in 0
£Y
£L/255 and 0
£X
£255. It means that a histogram of this data calculated one every byte will have its distribution cover a full histogram of underlying data, such a data with a wide spreading histogram will not ( empirically) guarantee a good achieve a data file with more compact histogram. For L different combinations, the X-Y coordinates of L points in a square of area bounded by
, ( as described above since L is much less than 65536 so that
must5 be considerably less than 255) and highly peaked occurrence concentrating on the first few integers. Also, it is not necessary to store this converting table because the inverse processor can reconstruct it from the stored 2D sorted histogram.
Examinations one some 2-D histograms of remote sensing images indicates that they are usually of the shape as shown in Fig.1. A break-point can usually be found at B which separates the curve in two parts. The first B most frequent combinations usually have a much compact shape with fast decreasing and generally cover a 60% to 80% of the area under the curve. It means that for two differential images 60% to 80% of the pixel combinations are of the first B combinations. Because we assign the (x,y) values or the code words symmetrically, that is the sequence of (x,y) would be (0,0), (0,1), ( 1,0) , (1,1), (0,2), (2,0),…(0. [
ÖL] + l}, ([
ÖL]+l,0), (1,2),(2,1(,… , and so forth, it is clear that if we first assign the (x,y) coordinates for the points in the square area bounded by 0
£X
£[
ÖB] + 1 and
0£Y£[ÖB]+1 and then the rest of points in the square of the area 0£X£[ÖL]+1 and 0£Y£[ÖL]+1 by the same rule , and the resultant data file must have a more higher peak on the first few integers, shown in Fig. 2. This is the actual processing in our system. And the processing is applied to every tow differential images generated by 2.1, and the break-point B is the index of the 2-D histogram array where its second derivative is minimized ( here the position of B is selected empirically based on the observations of many sorted histograms of combination values.
