# The Impact of KLT Coder on the Image Distortion in Astronomy

Authors: Petr  Pata
Journal: Advances in Astronomy
 (2015)

---

## 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.


---

## Body

## 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 2×2 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 ρ=1,…,M images with Gaussian probability density distribution: (1)xijρ,where i=1,…,N1 and j=1,…,N2 are column and row indices (i.e., pixel elements). These matrices are elements of vector space VN1N2 with defined scalar product(2)∑i=1N1∑j=1N2yijxij=Cxyover complex numbers. Then Cxy is a mutual energy of both elements. The symbol y is used for elements from dual vector space V~N1N2. The relation between VN1N2 and V~N1N2 is defined by relation (2). Orthogonality in defined vector space could be written as condition(3)∑i=1N1∑j=1N2yijxij=Cxyfor  orthogonal  elements0otherwise.When we can observe Cxy=1 for orthonormal elements, then we can find a set of independent elements(4)Φjlikj,l=1,k,l=1N1,N2which satisfy condition of the orthonormality(5)∑i=1N1∑j=1N2ΦjkilΦinjm=δlnδkmk,m=1,…,N1,l,n=1,…,N2.There, δln is well-known Kronecker delta. When the set is full, it can cover whole vector space VN1N2. Then the set is the orthonormal base of the vector space. The projection of the image to this orthonormal base is (6)ylk=∑i=1N1∑j=1N2Φkijlxjik=1,…,N1,l=1,…,N2.And inverse transform for the image reconstruction is(7)xmn=∑k=1N1∑l=1N2yklΦmknlm=1,…,N1,n=1,…,N2.The mean value could be obtained as(8)Exjiρρ=1M=∑ρ=1Mxjiρpdfxjim=1,…,N1,n=1,…,N2,where pdfxji is probability density function (pdf) of the distribution random image. The covariance matrix Ξ has elements(9)
Ξ
j
q
i
p
=
E
x
i
j
ρ
-
E
x
i
j
α
α
=
1
M
·
x
q
p
ρ
-
E
x
q
p
α
α
=
1
M
ρ
=
1
M
i
,
p
=
1
,
…
,
N
1
,
j
,
q
=
1
,
…
,
N
2
.When the pdfxji has Gaussian characteristic, the eigenmatrices of the covariance matrices are an optimal base in the vector space VN1N2 in terms of mean square error criterion. The characteristic equation of the covariance matrix is (10)∑i=1N1∑j=1N2ΞjnimΦjqip=βqpΦnqmpm,p=1,…,N1,n,q=1,…,N2.There, βqp are eigenvalues of the covariance matrix and are equal to energy contained in relevant spectral components. The set of base images Φjqipi,p=1;j,q=1N1;N2 can be sorted according to size of the eigenvalues βqpp=1;q=1N1;N2. 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 9×9 μm and 16-bit quantization depth.Figure 1
Size of eigenvalues versus their index for M7-300ff.dat image.Figure 2
Image from wide field camera M7 and Milky Way with many small objects (less than 10 pixels) and its autocorrelation function.


## 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 ofq=16 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 Figure2). 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 Figure3). These images are also mentioned as deep sky images (DSLI).The compressed image X⃡ with size Nx1×Nx2 can be decomposed by the operator R application into the set of M image submatrices with size N1×N2:(11)X⃡→Rxjiρi,j=1,α=1N1,N2,M,where, obviously, Nx1·Nx2=N1·N2M. 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 dimensionN1×N2 with full image covering.(ii)
Block with dimensionN1×N2 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 as(12)ξvu=∑p=1N1∑q=1N2xq+vp+uρ-Exq+vp+uαα=1M·xqpρ-Exqpαα=1Mu=-N1,…,N1,v=-N2,…,N2.The 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 VN1N2 can be found as N1=26–36 and N2=26–36 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 N1=N2=32 for used image classes.Figure 3
M42 nebulae with a satellite tray. Deep sky image from the BOOTES project (1024 × 1536 × 16 bits) with its autocorrelation function.Figure 4
Numerical complexity of the KLT, inverse KLT, and submatrices decomposition.


## 4. Image Coder Based on the KLT

Image compression coder can be composed from these parts (see Figures5 and 6):(i)
Application of theR 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 spaceVN1N2.(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.Figure 5
The integral transform part of KLT coder.Figure 6
The KLT coder for astronomical image compression.


## 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 ratio(13)CR=s0s1,where s0 is original image size s0=Nx1·Nx2·q and s1 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 withPSNR<25 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 as(14)MSE=1n·m∑i,j=1n,mPi,j-Qi,j2,where n·m is the total number of pixels of an image, i and j are indices of the image, and Pi,j and Qi,j are values of the pixels at location i,j of original image P and compressed image Q. The root mean square error (RMSE) is calculated as(15)RMSE=MSE.The 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 is(16)Br=Bmaxe-r2/2σ2+Bbgand the Moffat function is more general:(17)Br=Bmax⁡1+r/σβ+Bbg.There, r is a radial distance from the center of the stellar object (with maximal value Bmax), σ is a standard deviation, Bbg is a background value, B(r) is a pixel value at the radial distance r, 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 σx and σy (while considering the Gaussian function), where σx and σy are the distances from the center of the Gaussian at which the Gaussian is equal to e-0.5 of its central value. σx and σy 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 ofN1=N2=32 pixels each. The total number of realizations is M=1024. 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.Figure 7
Cumulative sum of the eigenvalues versus number of used KLT coefficients. Eigenvalues are sorted from the least to the most significant.Figure 8
The object position error (astrometry measurement) as a function of number of used KLT spectral coefficients.Figure 9
RMSE measurements of light image M7-300ff.dat (wide field camera) versus compression ratio.Figure 10
RMSE measurements of light image M17.dat (deep sky camera) versus compression ratio.Figure 11
Point spread function deformation (i.e., ellipticity) as a function of the compression ratio (wide field image).Figure 12
Point spread function deformation (i.e., ellipticity) as a function of the compression ratio (deep sky image).Astrometry analysis is also done for adaptive wavelet algorithm (JPEG2000) and for KLT coder. The comparison is shown in Figure13. 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.Figure 13
Comparison astrometry error for the KLT and JPEG 2000 coder.


## 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.



---

*Source: 484379-2015-11-24.xml*