Abstract
The data size of hyperspectral image is too large for storage and transmission, and it has become a bottleneck restricting its applications. So it is necessary to study a high efficiency compression method for hyperspectral image. Prediction encoding is easy to realize and has been studied widely in the hyperspectral image compression field. Fractal coding has the advantages of high compression ratio, resolution independence, and a fast decoding speed, but its application in the hyperspectral image compression field is not popular. In this paper, we propose a novel algorithm for hyperspectral image compression based on hybrid prediction and fractal. Intraband prediction is implemented to the first band and all the remaining bands are encoded by modified fractal coding algorithm. The proposed algorithm can effectively exploit the spectral correlation in hyperspectral image, since each range block is approximated by the domain block in the adjacent band, which is of the same size as the range block. Experimental results indicate that the proposed algorithm provides very promising performance at low bitrate. Compared to other algorithms, the encoding complexity is lower, the decoding quality has a great enhancement, and the PSNR can be increased by about 5 dB to 10 dB.
1. Introduction
Hyperspectral imaging technology is a major innovation in the course of development of remote sensing technology, which combines ground target spectra determined by material composition and the image reflecting ground target existing layout completely, and provides spectral information of dozens to hundreds of narrow bands while obtaining target spatial information for each pixel, thereby achieving a unity of image and spectra [1]. Hyperspectral remote sensing has been successfully applied in many areas such as geology, ecology, atmospheric study, and soil study and is playing an increasingly important role [2].
Hyperspectral image adds 3D spectral information based on 2D information about ground image. Hyperspectral data is a cube of a spectral image, whose spatial image dimension reflects the 2D spatial distribution of ground target features and whose spectral dimension reflects the spectral curve features. Spatial correlation exists in the spatial image dimension of hyperspectral image data and spectral correlation exists in the spectral dimension [3]. Accordingly, the difference between hyperspectral image compression and ordinary image compression is to remove both the spatial redundancies and the spectral redundancies. At present, the commonly used hyperspectral image compression methods include prediction-based, transformation-based, and vector quantization-based methods. For example, [4] used median prediction algorithm for hyperspectral intraband prediction and then used linear prediction and context prediction mixed algorithm for hyperspectral interband prediction, thereby implementing lossless compression of hyperspectral images. Penna et al. [5] used point extraction strategy to reduce the calculated amount of KLT and applied this improved KLT to the decorrelation of hyperspectral images, without having a significant impact on reconstructed image quality. Qian [6] proposed a fast vector quantization algorithm to improve the generation efficiency of codebook, which did not require full search, thereby greatly reducing the complexity.
In addition to the above conventional coding methods, fractal coding method characterized by its innovative ideas, high compression ratio, and resolution independence is considered as one of the most promising second-generation compression coding methods. Our previous studies of fractal video compression [7, 8] also demonstrate its good compression performance. Fractal coding techniques are based on the theory of iterated function systems (IFS) and are further developed by Barnsley in the early 1980s [9, 10]. The IFS theory is based on Banach’s contraction mapping principle, which states that a contractive transformation, defined on a complete metric space, possesses a unique fixed point or attractor. For the purpose of image compression, this idea translates into finding an optimal contractive transformation whose attractor closely approximates a given target image. This problem is widely known as the inverse problem in the fractal image coding literature. But this must involve manual intervention and it is extremely hard to find a contractive transformation with an attractor approximating the whole image. In the late 1980s, Jacquin developed a block-based fractal image compression scheme exploiting local self-similarities that are inherent in many real-world images [11] and made fractal image compression become automatic. In his scheme, the image is subdivided into a pair of simple and uniform partitions: a domain partition of larger subblocks, also known as parent subblocks, and a range partition of smaller subblocks, also known as child subblocks. A parent subblock is mapped into its corresponding child subblock using a geometric mapping, followed by a simple affine transformation, known as the gray-level map. Therefore, a digitized image can be stored as a collection of the position of the corresponding parent subblock and transformation coefficients for each child subblock. This generally requires much less memory, resulting in data compression.
We observe that there exist strong similarities between adjacent bands of hyperspectral images, so that we guess that fractal coding will be suitable for hyperspectral image compression. In this paper, we make a try to combine prediction and fractal coding for hyperspectral image compression. The algorithm first carries out intraband prediction for the first band which is treated as the basic band, and we extend the block-based fractal image compression scheme to two adjacent bands, which is to say each range block is approximated by the domain block in the adjacent band, which is of the same size as the range block to remove the spectral redundancies. The experimental results with AVIRIS (airborne visible infrared imaging spectrometer) hyperspectral remote sensing images indicate that the algorithm provides very promising performance at low bitrate.
The rest of the paper is organized as follows. The basis of fractal image coding is summarized in Section 2. The proposed hybrid prediction and fractal hyperspectral image compression algorithm is presented in Section 3. The experimental results are presented in Section 4. And finally the conclusions are outlined in Section 5.
2. Fractal Image Coding Basis
2.1. Fractal Image Coding Theory
The fundamental principle of fractal coding consists of the representation of an image by a contractive transform of which the fixed point is close to that image [12]. The image space is accepted as a complete metric space. By the contractive mapping theorem and collage theorem, the contractive transform is always possible within certain threshold. Original approach with IFS tries to find a number of affine mappings on the entire image, which is rather slow in terms of searching the contractive map function. Jacquin’s partitioned iteration function system (PIFS) [13] takes a different approach by finding the individual mappings for subsets of the images. Each PIFS contains a complete metric space and a series of contraction mapping as defined in this space. Fractal image coding just uses a PIFS to represent an original image so that the image after iterative decoding closely approximates the attractor of this PIFS as well as the original image. The contraction mapping coefficients constitute the fractal code of the original image.
The attractor theorem is as follows.
Attractor Theorem. Assuming is an IFS in the completely metric space , then we have the following.
Transformation defined by the following equation:is the contraction mapping of complete metric space , whose contractive factor iswhere represents the Hausdorff metric.
Compression transformation has a unique fixed point , so that it satisfies the following equation:
And the fixed point can be obtained by iteration, namely:where
Fixed point is called the attractor of IFS.
That is to say, for any set , iterative sequence converges to the attractor of PIFS. Its inverse problem takes the given original image as the attractor to solve its contraction mapping. It is very difficult to solve the inverse problem accurately, but collage theorem gives out a possible solution to this problem. That is, approximately minimize the distance between original image and attractor image by minimizing the distance between original image and its collage image .
The collage theorem is as follows.
Collage Theorem. Assuming is a complete metric space, given a set and a number , if we can find a , which makes thenwhere is Hausdorff metric and is attractor of this IFS.
2.2. Fractal Image Coding Implementation
The basic block-based fractal image coding algorithm flow chart is shown in Figure 1.
Firstly, divide the original input image into nonoverlapping blocks to form block set and overlapped blocks to form block set , with size of each block at and that of each block at . blocks and blocks are usually obtained by sliding window method.
Then, contract all blocks to the size of , that is, the same size as blocks, usually by downsampling or averaging the intensities of its four neighboring pixels,where represents the intensity of block at position .
Afterwards, contracted block is extended through 8 isometric transformations (including identical transformation; rotations by 90°, 180°, and 270° about the block center; symmetric reflections about vertical central axis, horizontal central axis, leading diagonal, and secondary diagonal). The extended domain pool is defined as .
In order to find the best matching block for each block, contraction mapping transformation is defined aswhere and represent spatial coordinates; represents pixel value; , , , and define one of the 8 isometric transformations; represents brightness adjustment factor with its absolute value smaller than 1; represents brightness offset factor.
Searching for matching block is a process of minimizing the following equation:where represents vector 2-norm.
Finally, use (9) and (10) to find the best matching block from the extended domain pool for each block and store the sequence number of isometric transformation, location of the best matching block, and corresponding brightness scale factor and offset factor . After fractal parameters of all blocks are stored, the total fractal encoding process is completed.
Fractal image decoding process is very simple, which can be completed by only inputting an arbitrary initial image (initial size of the image should be equal to that of its original image as required) and then iteratively applying the corresponding transformation to initial image based on the stored fractal parameters. Decoding process can be expressed as follows:where represents the th block of image after times’ iterations; represents the best matching block of and is derived from image after times’ iterations; is derived from the best matching block of the initial image; and represent the corresponding optimum brightness adjustment factor and offset factor, respectively; and represent spatial contraction transformation and isometric transformation, respectively. The symbol is the composition operator.
If reaches a predetermined number of iterations (generally less than 10), then stop iteration and output image . This is the decoded image, namely, approximation of the image to be encoded.
3. Hybrid Prediction and Fractal Hyperspectral Image Compression
Because each band of a hyperspectral image corresponds to the same location on the earth, hyperspectral data cube has a strong spectral correlation, whose spatial correlation is relatively weak. The basic fractal image coding algorithm does not use spectral correlation, which therefore cannot be directly used for compressing hyperspectral images. In order to take full advantage of spatial and spectral correlations of hyperspectral image, we propose a method for lossy compression of hyperspectral image based on hybrid prediction and fractal. The first band is intraband prediction encoded, and the remaining bands are interband fractal encoded. The difference is that each range block is approximated by the domain block in the adjacent decoded band, which is of the same size as the range block to exploit the spectral correlation. Then the prediction errors as well as fractal residual are further transformed, quantified, and entropy encoded, and the fractal parameters are also entropy encoded to improve the compression efficiency. The compression scheme divides each band into range blocks with the size of 16 × 16. The block diagram for this scheme is shown in Figure 2.
3.1. Intraband Prediction
Linear prediction has been studied in hyperspectral image compression [14–16]. In this section, we describe the intraband prediction design of our scheme. The main purpose is to remove spatial correlation as well as obtaining a high quality decoded reference band.
Hyperspectral images are usually texture abundant. For example, Figure 3(a) shows a simulated color IR view of an airborne hyperspectral data flight line over the “Washington DC Mall,” which has much detail. The sensor system used in this case measured pixel response in 210 bands in the 0.4 to 2.4 μm region of the visible and infrared spectrum. Bands in the 0.9 and 1.4 μm region where the atmosphere is opaque have been omitted from the dataset, leaving 191 bands. The dataset contains 1208 scan lines with 307 pixels in each scan line. The color image is made using bands 60, 27, and 17 for the red, green, and blue colors, respectively. Figure 3(b) is the first band of the hyperspectral image.
(a) The simulated color view with bands 60, 27, and 17
(b) The first band image
To predict more precisely, we design 9 linear predictors with different prediction directions, as shown in Figure 4. Eight prediction modes, as shown in Figure 4(a), are utilized for more suitably predicting textures with structures in the respective directions. The 16 samples of the 4 × 4 block, which is labeled as a–p, are predicted using the prior decoded neighboring pixels in the corresponding directions in the adjacent blocks labeled as A–Q, which is shown in Figure 4(b). In addition, the DC mode is used to predict smooth area, and each pixel in the shadowed block is predicted using the mean value of the neighboring pixels labeled as A–D and I–L, which is shown in Figure 4(c).
(a) 8 prediction directions
(b) Directional prediction
(c) DC prediction
The prediction mode with the minimum sum of squared difference (SSD) is chosen as the final prediction mode. The SSD is calculated as follows:where represents the source signal at position ; represents the reconstruction signal at position using the corresponding prediction mode.
3.2. Interband Fractal Coding
In order to remove the existing spectral redundancy in hyperspectral image, we have made some modifications to the basic block-based fractal image compression scheme. In our scheme, each range block is approximated by the domain block in the adjacent band, which is of the same size as the range block. And considering that the adjacent bands have similar spatial topology structures, as shown in Figure 5, thus isometric transformation process can be omitted. Statistics also show that the best matching block appears near block with a large probability [17], so it is not necessary to use full search. In our scheme, the search range of block matching is reduced to a small region centering the corresponding position of block, as shown in Figure 6. In addition, we adopt an adaptive block size instead of the original fixed size block partition method.
(a) The 25th band
(b) The 26th band
The specific steps for interband fractal coding are as follows.
Step 1. The current hyperspectral band image is divided into nonoverlapping blocks with the size of 16 × 16, and all blocks can cover the whole band image.
Step 2. Encode each block in turn.
The searching area is located in the adjacent decoded band, and with the same coordinate of the current block as the center, move 7 pixels (a variable parameter determined at the encoder) upward, downward, and towards left and right, to form the rectangular search area, as shown in Figure 6. Search the best matching block in the search area, with the MSE (mean square error) as the matching criterion:where represents the pixel value of the current block, represents the pixel value of the corresponding matching block, represents the number of pixels of the current block, is the scale factor, and is the offset factor. and are obtained as follows.
Solving the equation for ,
Setting to 0,
Similarly,
Setting to 0,
Substituting,
If the minimum MSE is less than T_16 (here is 10.5), the domain block can be considered as well-matched, and we can record the current IFS coefficients and continue to the next block; otherwise, go to .
The block is divided into two 16 × 8 subblocks. Search the best matching block for each 16 × 8 subblock, respectively, in the domain pool, and the size of searched matching block is 16 × 8. If the minimum MSEs of the two 16 × 8 subblocks are both smaller than T_16, the domain block can be considered as well-matched, and we can record the IFS coefficients of the two 16 × 8 subblocks and continue to the next block; otherwise, go to .
The block is divided into two 8 × 16 subblocks. Search the best matching block for each 8 × 16 subblock, respectively, in the search area, and the size of searched matching block is 8 × 16. If the minimum MSEs of the two 8 × 16 subblocks are both smaller than T_16, the domain block can be considered as well-matched, and we can record the IFS coefficients of the two 16 × 8 subblocks and continue to the next block; otherwise, go to .
This block is divided into four 8 × 8 subblocks for which, respectively, search for the best matching block also in the size of 8 × 8 in the search area, and if the minimum MSEs of all four 8 × 8 subblocks are smaller than T_8 (here is 8.0), we can record the IFS coefficients of the four 8 × 8 subblocks and continue to the next block; otherwise, we can further divide the 8 × 8 subblock which is not well matched into two 8 × 4 sub-subblocks or two 4 × 8 sub-subblocks or four 4 × 4 sub-subblocks just in the same way.
Step 3. Repeat Step 2 until all blocks of the current hyperspectral band have found the best matching blocks.
The flowchart is shown in Figure 7.
3.3. Residual Image and Fractal Parameters Encoding
In a general prediction encoding system, once the prediction is formed, it is subtracted from the original band, and the residual is usually further compressed using transformation encoding and entropy encoding. In our scheme, the residual image of the first hyperspectral band is obtained through subtracting the predicted band image by the original band image. In addition, the residual of fractal encoding is also processed in the same way, while it is ignored in the basic fractal encoding algorithm and most fractal literatures. The residual image of the fractal encoding is formed by
The residual image is DCT transformed, quantized, and Golomb entropy encoded. In parallel, the quantized data are rescaled and inverse transformed and added to the prediction band to reconstruct a coded version of the band which is stored for later interband fractal encoding. The fractal parameters are also Golomb entropy encoded.
4. Experimental Result
Two AVIRIS [18] hyperspectral data cubes, “Cuprite” and “Low Altitude,” derived from JPL are used. The size of the data cubes was originally 512 lines × 614 pixels with 224 spectral bands and 1087 lines × 614 pixels with 224 spectral bands, respectively. In this paper, a subset of them with 512 lines × 512 pixels with all 224 bands was tested. The experiments were performed on a computer with an Intel Core 2 Quad Q9300 @ 2.50 GHz CPU and 4 GB RAM. PSNR is used to measure the quality of the compression data. It is defined as follows:where Peaksignal is the maximum value in the data cube and MSE is calculated aswhere and are, respectively, the pixel values of the original and the reconstructed data of band at location , is the total number of lines in the data cube, is the total number of pixels per line, and is the total number of bands in the data cube.
Figure 8(a) shows the 17th band of “Cuprite” truncated from the original data cube for the experiment, and Figure 8(b) shows the decoded result image at 0.1 bpppb. As seen in Figure 8, there is almost no difference of subjective quality between the compressed image and the original image.
(a) Original image
(b) Decoded image
Figure 9 shows the PSNR comparison of “Low Altitude” encoded by the 3D-SPIHT [21] and the proposed algorithm. 3D-SPIHT coding is the three-dimensional extension version of the set partitioning in hierarchical trees (SPIHT) coding, which is based on the concentrative energy of wavelet coefficients, to exhibit symmetric branching in all three dimensions. More details can be found in [21]. It can be seen from Figure 9 that the PSNR of the first band by the algorithm of the article is higher than every other band, indicating that the prediction algorithm adopted obtains high encoding quality. The PSNRs of all bands of “Low Altitude” compressed by the proposed algorithm are higher than those in the 3D-SPIHT algorithm. It is calculated that the average value of PSNRs of all bands of “Low Altitude” compressed by the proposed algorithm is 60.07 dB, which is 15.91 dB higher than that of the 3D-SPIHT algorithm.
The PSNR comparison results at different bitrates are as shown in Tables 1 and 2. The AT-3DSPIHT [19] algorithm in Tables 1 and 2 is an improvement of 3D-SPIHT algorithm by means of asymmetric structure, having a longer tree along the wavelet axis than that of the traditional tree structure, and the PSNR performance of AT-3DSPIHT for hyperspectral image is better than 3D-SPIHT. More details can be found in [19]. Experimental results show that the PSNR of “Low Altitude” compressed by the proposed algorithm is slightly lower than that of AT-3DSPIHT at 0.1 bpppb, except that they are all higher than the contrast algorithms. In comprehensive consideration of all bitrates, compared with AT-3DSPIHT algorithm, the PSNR of the proposed algorithm is increased by 9.53 dB averagely, and compared with the algorithm in [20], it is averagely increased by 11.12 dB.
In [22], a compression algorithm of hyperspectral remote sensing image based on 3D wavelet transform and 3D fractal is proposed. The classical eight kinds of affine transformations in 2D fractal image compression are generalized to nineteen for the 3D fractal image compression. When the compression ratio on the Earth Observation-1 (EO-1) satellite hyperspectral image of China Dongting Lake area is 80 : 1, the PSNR is 36.78 dB, while the PSNR of our proposed scheme at the same compression ratio is 41.78 dB. We can achieve a compression quality improvement of 5 dB.
To evaluate the encoding complexity of the proposed algorithm, we compare the encoding time of the proposed algorithm with that of Lucana Santos’ [23]. The corresponding encoding time by the two schemes for “Yellowstone_sc0” and “Yellowstone_sc3” is recorded in Figure 10, which is displayed in seconds/band. Both scenes have the size of 512 lines × 680 pixels with 224 bands. From these results, we can see that the proposed algorithm has lower complexity and our scheme saves total 94.1% encoding time on average.
(a) Yellowstone_sc0
(b) Yellowstone_sc3
5. Conclusion
A lossy compression algorithm of hyperspectral image based on hybrid prediction and fractal is proposed. Intraband predictors with different directions are designed to remove the spatial correlation of the first band, and a high decoding quality can be obtained, which provides a good basic band for the following fractal encoding. The remaining bands are encoded by interband fractal encoding. The spectral redundancy can be removed since each range block is approximated by the domain block in the adjacent band, which is of the same size as the range block. Adaptive block partitions as well as a local search are also adopted in the fractal encoding. Finally, the residuals of prediction and fractal encoding are DCT transformed and Golomb entropy encoded. The fractal parameters are also Golomb entropy encoded.
Experimental results show that the algorithm achieves effective compression of hyperspectral image at low bitrate. The decoded quality is greatly improved, and the PSNR performance is increased by about 5 to 10 dB. Besides, the encoding complexity also has a great advantage. The proposed hybrid prediction and fractal hyperspectral image compression algorithm is a breakthrough to the traditional algorithm framework of hyperspectral image compression, providing a commendable solution for lossy compression of hyperspectral image.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This project is funded by the National Natural Science Foundation of China (NSFC) under Grants nos. 61375025, 61075011, and 60675018 and also by the Scientific Research Foundation for the Returned Overseas Chinese Scholars from the State Education Ministry of China. The authors would also like to express their appreciations to the reviewers for their thorough review and very helpful comments, which help in improving this paper.