Abstract
Presented paper is devoted to the application of Karhunen-Loève transform (KLT) for compression and to study of KLT impact on the image distortion in astronomy. This transform is an optimal fit for images with Gaussian probability density function in order to minimize the root mean square error (RMSE). The main part of the encoder is proposed in relation to statistical image properties. Selected astronomical image processing algorithms are used for the encoder testing. The astrometry and point spread function distortion are selected as the most important criteria. The results are compared with JPEG2000 standard. The KLT encoder provides better results from the RMSE point of view. These results are promising and show the novel approach to the design of lossy image compression algorithms and also suitability for algorithms of image data structuring for retrieving, transfer, and distribution.
1. Introduction
Novel robotic imaging systems acquire huge amount of data during their service [1]. This data is necessary for archiving and distribution among servers and users. Therefore, searching for suitable compression standard is still rather important task [2, 3]. There are many compression standards matched to the human visual system (HVS) for multimedia purposes [4]. These methods are optimized for subjective image quality. The most popular standards are JPEG2000 (based on the wavelet filtering) and classical JPEG with discrete cosine transforms. Multimedia compression algorithms have different quality criteria. The HVS is fundamental for these standards [5]. In particular, they are human visual perception properties such as spatial and temporal resolution, sensitivity to brightness difference, and spectral response. Astronomical images are devoted to completely different purpose and processing tasks. Therefore, this novel standard may not be suitable for an application aimed at human observer. This is typical for astronomical image compression in astronomy, medicine, and other scientific applications [6].
The lossless standards are often preferred for compression in astronomy. The popular algorithm is FITSPRESS developed at the Center for Astrophysics, Harvard. The FITSPRESS algorithm uses Daubechies-4 filters of the wavelet transform and application of the run length encoding and Huffman code [7]. The second standard HCOMPRESS has been developed at the Space Telescope Science Institute (STScI, Baltimore). The HCOMPRESS is used for distributing archive images from Digital Sky Survey DSS1 and DSS2. The Haar transform with blocks of pixels is used for this standard and it is extremely fast [8]. The multiscale approach is also suitable for lossless image compression. This principle is included in the astronomical context coder (ACC) [9]. Lossless algorithms have limited capability of achievable compression ratio (up to 5 : 1) according to noise level in image data and limited redundancy.
The lossy algorithms offer higher compression ratios than lossless approach. The removal of irrelevancy from image can bring an irreversible distortion of compressed image. It has a direct influence on distortion of data and therefore it is necessary to study the distortion nature. The PHOTZIP is a lossy method based on the image function modeling. The compression ratio can be driven by the setting of acceptable astronomical measurement error [10]. Popular modern image compression standards (JPEG2000, JPEG-LS, or JPEGXR) offer also high compression ratios for wide range of image classes including astronomical applications [11, p. 200] [12]. These standards are designed with very good performance for common images [13]. They can take advantage from special properties of selected special image types [14]. These classes include biomedicine, security, and satellite applications. The special effort has been spent for robotic telescope and other systems with large amount of acquired image data.
The Karhunen-Loève transform (KLT) is the integral transform with optimal data reduction properties. The KLT offers the best data fit with root mean square error minimization for image data with Gaussian probability distribution. The KLT has also close relation to PCA or Hotelling transform used in image processing and pattern recognition [15]. The main advantage of this transform is a signal decomposition into optimal decorrelated components in KLT vector space. The KLT could be extended for applications in image processing [16] and components decorrelation could be used in image compression. The optimal KLT bases are signal dependent and it is necessary to construct base system for each image separately for the optimal Gaussian decorrelation. These characteristics give high potential of the KLT for signal and image processing based on the decorrelation. Compression is one of these tasks. The optimal decorrelation is highly efficient for hyperspectral imaging especially [17]. The lossy compression algorithms can be used for large image archives and intelligent management of image data [3, 18].
2. Karhunen-Loève Transform
The Karhunen-Loève transform (KLT) is well-known integral transform for Gaussian signals processing [15]. Let us assume that there is a set of images with Gaussian probability density distribution: where and are column and row indices (i.e., pixel elements). These matrices are elements of vector space with defined scalar productover complex numbers. Then is a mutual energy of both elements. The symbol is used for elements from dual vector space . The relation between and is defined by relation (2). Orthogonality in defined vector space could be written as conditionWhen we can observe for orthonormal elements, then we can find a set of independent elementswhich satisfy condition of the orthonormalityThere, is well-known Kronecker delta. When the set is full, it can cover whole vector space . Then the set is the orthonormal base of the vector space. The projection of the image to this orthonormal base is And inverse transform for the image reconstruction isThe mean value could be obtained aswhere is probability density function (pdf) of the distribution random image. The covariance matrix has elementsWhen the has Gaussian characteristic, the eigenmatrices of the covariance matrices are an optimal base in the vector space in terms of mean square error criterion. The characteristic equation of the covariance matrix is There, are eigenvalues of the covariance matrix and are equal to energy contained in relevant spectral components. The set of base images can be sorted according to size of the eigenvalues . The size of eigenvalues is shown in Figure 1 for a wide field image example. M7-300ff.dat is the image from the BOOTES system (see Figure 2). These images have size 1536 × 1024 pixels with pixel size of μm and 16-bit quantization depth.
3. Image Data
The proposed approach is designed for images obtained from automatic systems BOOTES [18], MAIA [19], WILLIAM [20], and GLORIA [1]. BOOTES (Burst Observer Optical Transient Exploring System) is the first Spanish robotic telescope and it is a unique system with ambitious target continuous sky observation. The BOOTES has four main stations (Huelva and Málaga in Spain, Blenheim in New Zealand, and Kunming in China). The first light was obtained in 1998 and main aim is focused on the observation of the extragalactic objects and GRB optical transients.
The MAIA (Meteor Automatic Imager and Analyzer) is a video system for observation of weak meteors [21]. The two MAIA stations are equipped with high quality Ethernet cameras with image intensifiers with nonlinear response function. The special MAIA software is designed for automatic processing video stream. The WILLIAM (wide field all-sky image analyzing monitor) system is a low cost UWFC (ultrawide field camera) experiment. The main purpose of the WILLIAM system is autonomous weather monitor for robotic telescopes. The GLORIA (Global Robotic-Telescope Intelligent Array) is free and open-access network of robotic telescopes; for more details see [22].
The proposed image coding scheme has been improved on the database of more than 200 images from BOOTES, MAIA, and WILLIAM. Images have precision of bits and are in nonbinning mode. Dark frames have not been extracted from original images. RAW image artefacts especially hot pixels and cosmic trays have very narrow response of the autocorrelation function and they are compressed with difficulty. The images in database have been divided into two main classes according their response autocorrelation function:(i)Wide field images (images from the short focal length cameras): the autocorrelation function is very narrow (see Figure 2). Typically these images are from WF and UWF cameras. The object FWHM is few pixels only.(ii)Deep field images (images obtained from the primary focus of the telescope): these systems have a high spatial resolution. The FWHM of the autocorrelation function is more than 10 pixels (see Figure 3). These images are also mentioned as deep sky images (DSLI).The compressed image with size can be decomposed by the operator application into the set of image submatrices with size :where, obviously, Let these submatrices be assumed as realizations of random process in Hilbert vector space [11]. The decomposition operator can be realized in these principal modes:(i)Block operator similar to classical JPEG principle with block dimension with full image covering.(ii)Block with dimension centered in detected objects.(iii)Wavelet band decomposition.(iv)Decomposition operator related to irregular areas around detected objects in the image.The interpixel correlation is an important criterion for block size determination. The 2D autocorrelation for the block number can be calculated asThe 2D autocorrelation curves of typical images from archive are shown in Figures 2 and 3. The block size is a compromise among coder efficiency, computational complexity, and image data decorrelation. The optimal dimension of the space can be found as and pixels for images with size 1536 × 1024. This result has been determined from the autocorrelation function shape and numerical cost measuring (see Figure 4). The numerical cost can be divided into two parts. The first part is time for image decomposition according to (11) and the second one is numerical complexity for KLT and inverse KLT calculation based on (6) and (7). The relation between time cost and size of the decomposition submatrix is shown in Figure 4. The optimal size has been chosen equal to for used image classes.
4. Image Coder Based on the KLT
Image compression coder can be composed from these parts (see Figures 5 and 6):(i)Application of the decomposition operator.(ii)KLT spectrum calculation.(iii)Quantization of spectral coefficients.(iv)Lossless coding for redundancy removing.A special order of the spectral coefficients has been derived for easy image reconstruction. The data amount reduction has been performed for the designed coder in these steps:(i)Selection of the set basis vectors of the proposed vector space .(ii)Decorrelation of image data in the KLT spectrum.(iii)Reduction of the spectral coefficients number.(iv)Nonlinear quantization of spectral coefficients.(v)Organization of the data stream for the lossless compression (Huffman scheme).(vi)Calculation of the different matrix for image error estimation.
5. Image Distortion Analysis
The KLT sorts spectral components according to their importance with the mean square error minimization. This arrangement offers minimal image distortion and with sophisticated data stream organization has the minimal influence on the astrometry and photometry measurements. These criteria have been chosen for the image distortion analysis. The KLT decomposition has been found as an effective approach to astronomical image compression. The compression efficiency can be evaluated by ratiowhere is original image size and is compressed image size, respectively.
The most commonly used method of the distortion evaluation is also well-known PSNR (Peak Signal to Noise Ratio). The PSNR measures image distortion, usually quoted in decibels, in a logarithmic scale. The criterion has a limited approximate relationship with the perceived errors noticed by HVS (human visual system). The compressed image with PSNR higher than 30 dB has defects indistinguishable from the original image. The image with dB is not often acceptable for subjective image quality point of view. This approach is more related to subjective methods of quality measurement in the multimedia applications.
Besides the PSNR, there are MSE and RMSE criteria; they are often used to determine the quality of compressed images. The mean square error (MSE) can be expressed aswhere is the total number of pixels of an image, and are indices of the image, and and are values of the pixels at location of original image and compressed image . The root mean square error (RMSE) is calculated asThe second group of the image distortion metrics includes methods based on data-evaluation algorithms. These methods are based on stellar profile fitting (i.e., point spread function, PSF, and defects), aperture photometry, astrometry (position error), or successful detection of stars [23–25]. There are two common functions for fitting star profiles: Gaussian and Moffat. The effort is to match a star profile as close as possible to the Gaussian or Moffat profile. The Gaussian function isand the Moffat function is more general:There, is a radial distance from the center of the stellar object (with maximal value ), is a standard deviation, is a background value, is a pixel value at the radial distance , and is a Moffat coefficient. Center of a star usually has a profile closer to the Gaussian function. Peripheral parts are closer to the Moffat profile. So the ideal fitting function combines Gaussian and Moffat. The ultrawide field systems have a point spread function different from the ideal Moffat fit. The space variant approach is necessary to use for such systems [20]. The point spread function of these systems has different shape from the Moffat and Gauss approximation. The stellar profile has to be fitted by sophisticated Zernike polynomials which are more general for systems with field dependent aperture. The shape of objects can be described as an ellipticity parameter. The ellipticity is defined by and (while considering the Gaussian function), where and are the distances from the center of the Gaussian at which the Gaussian is equal to of its central value. and are the major and minor half-axis of the ellipse, perpendicular to each other (see Figures 11 and 12).
6. Results
The set of test images was chosen from the DEIMOS database [26] and included more than 200 images from BOOTES, WILLIAM, and MAIA experiments. There were three typical classes of images in the set. The first group contains deep sky images with high angular resolution (BOOTES) and the second ultrawide field (WILLIAM) and special sequences from system with nonlinear response function (MAIA).
Each image has been divided into blocks with size of pixels each. The total number of realizations is . The covariance matrix has been computed using abovementioned set of equations. The eigenvalues of covariance matrix have been used as a good estimation of the significance of relevant spectral coefficients. Eigenvalues express energy proportion in the corresponding spectral coefficients as it is shown in Figure 7. The RMSE of the object position after less significant spectral coefficients omitting is shown in Figure 8. More than 80% of the KLT spectral components are redundant from RMSE point of view. Their omitting leads to object position error less than 0.5 pixels. The adaptive postprocessing techniques can be used for the data and accuracy of astrometry measurement enhancement. The RMSE distortions are shown in Figures 9 and 10 for wide field and deep sky images. There, RMSE in levels is plotted (16-bit image) as a function of the achieved compression ratio. As one can expect, RMSE is a significant criterion for astronomical image compression. This is a different result from multimedia compression approaches.
Astrometry analysis is also done for adaptive wavelet algorithm (JPEG2000) and for KLT coder. The comparison is shown in Figure 13. Application of the KLT coder brings smaller influence on the object position than JPEG2000 for the same compression ratio. It can be explained as an effect of the wavelet filtering. This method smooths the signal with characteristic edge blurring. Wavelet effect distorts the point spread function of the system and shape of objects. Classical object fits based on Moffat function are not sufficient. The KLT has distortion with RMSE minimization and it has smaller influence on object shape as is shown in Figures 11 and 12.
7. Conclusion
The proposed coder based on the KLT can be found as a good alternative for known compression algorithms (JPEG and JPEG2000). The KLT coder has extensive computational requirements for eigenvectors calculation of covariance matrix. It can be improved using the KLT coder in modified arrangement. This approach is based on the calculation of eigenimages for typical image classes optimized for used imaging system. The KLT coder has smaller influence on object shape (PSF) than frequently used algorithms based on wavelet transform. Further improvement of technique could be achieved by sophisticated filtering methods and proper image data organization.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work has been supported by Grant no. 14-25251S “Nonlinear imaging systems with spatially variant point spread function” of the Czech Science Foundation.