Influence of lossy compression techniques on processing precision of astronomical images

2005 IEEE International Symposium on Signal Processing and Information Technology

Influence of Lossy Compression Techniques on Processing Precision of Astronomical Images Pavel Hanzlík, Petr Páta, Jaromír Schindler, Stanislav Vítek Department of Radioelectronics, Czech Technical University of Prague, Prague, Czech Republic Abstract - Compression of image data from Burst Observer and Optical Transient Exploring System (BOOTES) is discussed. BOOTES project is an international robotic telescope for optical transient of Gamma Ray Bursts (GRB). The statistical distribution of image functions in astronomical images from wide field and deep sky cameras is compared with Gaussian and Laplacean probability density function (pdf). The comparison of two irrelevancy reduction methods is presented from a scientific (astrometry and photometry) point of view. First one is based on a statistical approach to data compression and it is suggested from the Karhunen-Loève transform (KLT) with uniform quantization in spectral domain. Second technique is derived from wavelet decomposition with adaptive choosing of used mother wavelet. Keywords – astronomical image data compression, Karhunen – Loève expansion, adaptive wavelet decomposition, BOOTES project

I. INTRODUCTION There is a rapid development of new scientific image detection systems based on CCD or CMOS technology during last 10 years. The new high precision CCD sensors allow new applications in the geoinformatic (remote sensing of Earth), astronomy and satellite branches. Older photographic plates are substituted by CCDs with high spatial resolution and quantization depth. Therefore the costs for data transfer, processing and archiving are increasing. The MegaCam camera of Canada-French-Hawaii telescope (CFHT) can work with spatial resolution up to 16k x 16k pixel in 32 bits precision [4] or robotic telescopes in full automatic mode. There is no exception having an image in size of higher than 2 MB. The year data amount can be estimated as 632.8 GB in full mode for BOOTES robotic system [2]. Considering the limited capacity of storage media, an efficient data compression algorithm has to be applied. Loosless compression algorithms are often used in scientific applications but their efficiency is limited. Maximal achieved compression ratio depends first of all on the data type and amount of image signal entropy. The most efficient loos-less algorithms are the Huffman and arithmetic coding. Typical compression rates of loos-less algorithms are from 1:1.5 to 1:5 [1].

0-7803-9314-7/05/$20.00©2005 IEEE

Fig. 1. The BOOTES-1A, B station in the southern Spain (Mazagon).

The second possibility is usage of more efficient lossy compression techniques characterized by decorrelated parameters. Typical examples of this possibility are JPEG or JPEG2000 standards but the data impairment has to be taken into account. It is necessary to consider whether the algorithms optimized for multimedia applications and human vision are suitable for scientific image data compression. JPEG standard is recognized as a not suitable one because its quantization tables are defined for color image signal with maximal depth 12 bits per each color channel, while scientific (e. g. BOOTES) system produces 14 - 16 bit gray scale images. Application of special compression algorithms based on special properties of wavelet, fractal or Karhunen Loève expansion (KLE) [11] seems to be a better solution. Comparison of Karhunen - Loève expansion and wavelet approach for compression of image data acquired during BOOTES experiment is evaluated in this study. 2. BOOTES PROJECT The BOOTES project (Burst Observer and Optical Transient Exploring System) [5] is a system for the sky monitoring and searching the optical transient of GRB (Gamma Ray Bursts) without a human interaction. The system consists of two parts, BOOTES - 1 and BOOTES – 2. First one is located in the south Spain (Mazagon near the Huelva city) and has been in full operation since July 1998 (see Fig. 1). The second station is in eastern part of south Spain (La Mayora). The

346

3. IMAGE DATA PROCESSING Image data that can be acquired by watching the space with astronomical instruments are often distorted by a variety of factors. The main factors are the atmospheric distortion (scintillations), optical system inaccuracy (aberration, diffraction), distortion by electronic components (primarily image sensors) and at last but not least also image compression methods. The computational technology rise in last few years allows us to minimize the technology impacts of any inaccuracy on one side, but also, on the other hand, the resolution and data volumes are increasing rapidly. That on the contrary brings more problems with storing and post processing of data. As an example, let’s introduce the all-sky BOOTES monitor. Every minute it brings a frame in size of 4096x4096 pix. After the loosless compression (gzip) it introduces 6.5 GB of image data per night. Moreover its prepared an upgrade of BOOTES1A telescope for observing of Swift error-box - two cameras mounted on telescope will produce about 50 GB per night in the case of nice weather.

Fig. 2. Input image data from BOOTES project., 1024x1536x16 bits a) (above) Image from wide field camera M7 and Milky Way with many small objects (less than 10 pixels). b) (bellow) M42 nebulae with a satellite tray. Image has been obtained from DEEP SKY camera of BOOTES.

whole system was completed in July 2001. The BOOTES was the first Spanish robotic telescope [3]. The large base distance (240 km) between the stations is important for false phenomena recognition. The phenomena comprise images of near-Earth objects with a distance shorter than one million km, reflections from the satellites surface, meteors, cosmic rays effect, etc. The wide field cameras of the system has field of view 16x11 degrees with an observation depth up to 14 mag. The BOOTES stations are equipped with Paramount, Losmandy and Meade robotic telescopes, SBIG and Apogee CCD cameras and a unique spectrograph. The most important benefit of the project is a possibility of optical transient of GRB detection. This observation is in a close collaboration with other specialized space instruments – for example INTEGRAL, HETE and SWIFT satellites, etc. Discovery of new long variable stars and active galactic nuclei observation can be also considered to be as important outcomes of BOOTES work.

As said in the former parts, main goal of this article is a study of impacts of lossy image compressions. While still remaining in the branch of astronomy, we may use several different approaches – the best may be the evaluation of exactness and efficiency of algorithms used in astronomy to detect various types of objects, divide them to appropriate type and their evaluation and qualification. It means namely that astrometry deals with the positions of stars and other (stellar) objects (see Fig. 2), their distances and movements, and photometry, the process of determining the apparent magnitude (relative brightness) of stellar object. Photometric methods can be divided into the two main groups: methods deal with Point Spread Function (PSF) fitting and (CCD) aperture photometry. The brightness profile of the star in an image (PSF in fact) is Gaussian-wise. There are basically two functions, which are used in astronomy to model the star. The first one is a Gaussian function B(r ) = Bmax e

−r 2 R2

+ Bbg

(1)

where Bmax stands for maximal brightness – at the radial distance r = 0 , Bbg is a background brightness and R stands for characteristic radius. A more exact estimation of the object profile can be done by using the Moffat function [13] B(r ) =

Cmax r 1+   R

β

+ Bg

(2)

which by appropriate validation of the coefficient β allows us a very effective fit of the object profile. Since this type of image is very specific and basically the only object carrying he energy in higher spatial frequency (star) is well defined by Gaussian function, the adaptive estimation respecting the star model should be highly effective. To make the fit even more

347

efficient, it should have been in ideal case a n-order polymonic function. But the fit realization would be highly dependent to quality of reconstruction of the image data. We are using an approximation by combining the Gaussian impulse and Moffat function. 4. DERIVATE STATISTICS – GAUSSIAN APPROACH OPTIMALITY Fig. 3. Histogram for marginal pdf in a log plot, natural image (left) and tested scientific image (see Fig. 2 a).

We have studied the BOOTES 2 image data in order to evaluate their probability density distribution function and align the facts from the evaluation with a compression method which would be suitable for the project. 4.1 SINGLE PIXEL STATISTICS

In accordance with Huang [8] we show that being i, j the coordinates of the image and I(i, j) the luminosity function of the pixel, the probability density function (pdf) can simply be defined as a histogram of f(i,j), where

∑ (I f (i, j) = I (i, j) -

(i, j) )

i, j

i max j max

,

(3)

what only means we also take in mind getting rid of the mean value of the image intensity. In an optimal case, in astronomical image compression, we are looking for Gaussian probability density function, because in accordance with the usual notation Ν(m, σ2) we may model the probability density function of a Gaussian variable with variance σ such as f W (k ) =

1

σ 2π

e −k

2

/( 2σ 2 )

.

D = I (i, j) - I (i, j+1)

Obviously the difference between a natural image and an astronomical image shows small, but in terms of optimality, important detail (see Fig. 3). These marginal densities are rather similar for a defined group of images, with variance most distorted around the graph tails. We can watch a similar character at a defined set of images, smoother character for high resolution images i.e. astronomical images. The Gaussian character (as viewed from the relation (4)), is a turned parabolic shape in linear Cartesian coordinates. Changing to logarithmic, we get an optimal curve (see Fig. 4). It can be shown that since the character of the image data would nearer to Gaussian curve, the KLT will decorrelate the luminosity probability function more efficiently [6]. We can define two models, which align very well with the distinct peak at zero derivation and concave tails Model 1: f ( x) =

1 1 Z (1 + x 2 / s s ) t

1 − x/s α e Z For these models, Z is fixed since the integral both models then have two free parameters. known as the t-distribution and model 2 is Laplacean distribution where parameters s, α dependent on variance and kurtosis

Model 2: f ( x) =

(4)

Truly to say, in nature there won't be many images showing even partially Gaussian or any other basically known distribution character that's why it is very uneasy to find a model for single pixel statistics [7]. In astronomical images data, we tested a set of images in order to evaluate their common statistical character, what may help us to use an optimal compression method. In case of a Gaussian probability density function, Karhunen-Loève transform (PCA) would be the optimal transform. 4.2 MARGINAL DISTRIBUTIONS The probability density function (pdf) can also be formed through a histogram of adjacent pixel derivations (or differences). This directly connotes a dependence on orientation and need of defining the span of that adjacency i.e. correlation distance. The mean of the image intensity is also subtracted in this case.

(6)

κ= σ2 =

Γ(1 / α )Γ(5 / α ) Γ 2 (3 / α ) s 2 Γ(3 / α ) . Γ(1 / α )

(7) (8) of f(x) = 1, Model 1 is generalized are directly

(9) (10)

While having X a random vector on ℜ , we can define with µ (mean) and σ2 (variance) of X the kurtosis (κ) κ=

Ε( X − µ ) 4

σ4

and skewness (Ѕ) as follows

348

(11)

Fig. 4. a) – left Gaussian function in logarithmic coordinates b) - right Astro image form Fig. 2a) distribution approximated by Model 2.

S=

E( X − µ ) 3 .

(12)

σ3

Fig. 6. Block diagram of an adaptive wavelet encoder.

curve of an image matrix from BOOTES archive is shown in Fig. 5. Relation between bases in vector can be expressed | y lk ] =

N1 N 2

∑∑ | Φ

jk il

k = 1, L , N 1

] | x ij ]

l = 1, L , N 2 .

i =1 j =1

(14)

Equations for inverse transform can be derived applying orthogonality properties of the base | x mn ] =

∑∑ [Φ

The Karhunen-Loève expansion is an optimal transform of statistical data from mean square error point of view. The base vectors are composed from eigenvectors of covariance matrix of evaluated images. The important presumption for our analyze is a statistical correlation among image pixels. It can be considered as satisfied for a special image data – for example taken with the same spatial resolution, optics and similar shape of image function. We consider BOOTES image data to be sufficient. t Let X is an input image matrix with dimension N x1 × N x 2 .

Rt is an operator conducting decomposition of image matrix X into M sets of image submatrices with dimension N1 × N 2 α

N1 , N 2 , M

(|

(15)

x nm

M

i j

m n

=1

)

M M =1 )} ρ =1

(16)

where m, n are indexes within the entire base matrix and r , s describes the entire base matrix. 6. WAVELET APPROACH Image compression based on wavelet transform is nowadays very popular because of its nice properties. Discrete dyadic wavelet transformation is given by the equation

[ ] ∑ f [m]ψ

Wf n,2 j =

N −1

* j

[m − n]

(17)

m =0

where ψ

where, obviously, N1 .N 2 = N .N .M . The derived submatrices can be assumed as N1 × N 2 dimension realizations of random process in Hilbert vector space [11]. The 2D autocorrelation

ψ j [n] =

,

{α [ x |}α ] ρ − E {| x ]α }α

Ξ injm = E{( ρ [ x ij | − E

(13)

i , j =1,α =1

n = 1, L , N1 m = 1, L , N 2 .

|| y lk ]

Elements of the covariance matrix of the random process [11] can be determined as

5. KARHUNEN-LOÈVE EXPANSION

] 

mk nl

k =1 l =1

Fig. 5. The 2D autocorrelation curve of 1M11.00.dat.

t R  i X → |xj 

N1 N 2

349

j

is the wavelet on the resolution level j ∈ Z .

 n  Ψ j  2  2 1

j

(18)

Fig. 7. Original image of the star (left), Image with compression ratio 1:150, quantized to minimize the MSE (right).

and Ψ is so called mother wavelet. In real application we choose j ∈ {0,1,2,...i} and all remaining resolution levels are aggregated with

[

] ∑ f [m]φ N −1

Lf n,2 i +1 =

* i +1

[m − n]

Fig. 8. Comparison of impact of irrelevancy reduction for Adaptive wavelet algorithm (a – above), DCT (JPEG) (b - left), and DKLT (c - right) based coder.

(19)

m =0

where

φi +1

φ ï +1 [n] =

( [ ]

is a scaling function for resolution i + 1 1 2

i +1

Q j ,d ,b Wf n,2 j

 n  Φ i +1  2 

(20)

and Φ is so called father wavelet. Father wavelet is explicitly given by the mother wavelet and vice versa [10]. It is possible to construct such a father wavelet that the inverse wavelet transformation of Wf n,2 j j∈{0,1...i } and Lf n,2 i +1

[ ]

[

]

will describe original input function f [n] without the loss of any information. The efficiency of such a compression is given by selection of suitable mother wavelet. Than the wavelet transform represents the input image with the few non-zero coefficients. However, almost every image (not only astronomical) is highly non-stationary signal and a single wavelet cannot optimally correlate with its whole area. This leads to an algorithm, where the wavelet is determined according to very local properties of input signal. Orthogonality of the basis is not crucial so this freedom leads to combinatorial explosion. Since the compressed output stream cannot be overloaded by the information about wavelet change, the causality for the decoder must be added instead. Lossless compression ratio varies according to type of input image and the basis estimator effectiveness approximately from 3 to 5. This ratio is in most cases unsatisfaction so the quantization of wavelet coefficients is involved [12]. The quantization step makes more difficult to preserve the causality for the decoder thus the basis estimation process is less effective. The approximation f ' [m] of an image f [m] after quantization will be build from new set of coefficients

j∈{0,1...i

)

( [

i +1 } and Q j ,b Lf n,2

])

(21)

where Q j ,d ,b is a quantization operator which is usually different for each resolution level j , direction d in which the DWT is applied to an input image and for the block of coefficients defined within one resolution level and the DWT direction [9]. The image is then reconstructed with some loss of quality. It is clear that there is a lot of possible ways to quantize the coefficient in order to obtain desired compression ratio. But the only one setting is optimal with respect to distortion of scientific information included in the input image. 7. RESULTS AND MEASUREMENTS Images M7-300FF.st8.dat and M42_d03.dat (Fig. 2) acquired from BOOTES-1A system has been chosen as a testing image signal. The image has been taken with a short and narrow focus. Compensation of dark frame and flat field correction have been carried out. According to the very short focal length the size of the objects (stars) do not exceed 10 pixels (see Fig. 2a). Therefore the quality of applied compression algorithm is very important. The comparison measurement was made for DCT (JPEG), suboptimal KLT and adaptive wavelet. Distortion of image function has been impressed by irrelevancy reduction and comparison of these defects has bee shown in Fig. 8. The discrete cosine transform (JPEG) can be excluded as unacceptable according to characteristic blocking effect and removing of faint objects. Suboptimal Karhunen–Loève Expansion (KLE) is also the integral transform in discrete form indeed with more efficient spectral distribution. The characteristic block domain is noticeable without distortion of

350

the faint object profile for a linear quantization function in spectral domain (Fig. 8c). Higher compression ratios of adaptive wavelet transform required coarse quantization also in lower resolutions of wavelet coefficients, where the support of the wavelets was comparable to the objects sizes. This leads to objects distortions and artifacts characteristic to wavelet approach to data irrelevancy reduction.

3,5 Suboptimal KLT

Position error [pixel]

2,5

ACKNOWLEDGEMENT This work has been conducted at the department of radioelectronics of the Faculty of Electrical Enegineering of the Czech Technical University in Prague in the frame of the research project "Qualitative Aspects of Image Compression Methods in Multimedia Systems" and has been supported by grant No. 102/02/0133 of the Grant Agency of Czech Republic.

1

0,5 0 0

50

100 Compression ratio

150

200

Fig. 9. Error of the astrometry position measurement for the suboptimal Karhunen - Loève expansion and adaptive wavelet transform.

8. CONCLUSION

There are many possibilities how to quantize the wavelet coefficients and it would be possible to implement an adaptive quantization if we are able to detect objects we are focused on. This would improve the compression ration while maintain the scientific information of an image. KarhunenLoève expansion of image data offers an alternative compression method to the algorithms based on fractal and wavelet transforms. Disadvantages of this method are extensive computational requirements coming from necessity of solving of characteristic equation in the coder. Further improvement of the technique could be achieved by sophisticated filtering methods and suitable image data organization.

2

1,5

The object position measurement has been chosen as criterion for comparison of Karhunen -Loève expansion and adaptive wavelet transform. The astrometry has been done by fit of the Moffat function (2) and results are shown in Fig. 9. The acceptable deviations of object position are smaller than 0.5 pixels and compression ratio up to 50:1 for adaptive wavelet transform and to 100:1 for suboptimal Karhunen-Loève transform.

Since the astronomical image data nearer to Gaussian character, the Karhunen-Loève transform would be the most appropriate compression method for the data compression. The two defined models have a capability to either fit into the Gaussian character or conform variably to any differences from the optimum. If the character would be strictly Gaussian, we would be able to decompose the image into few decorrelated vectors with no non-diagonal elements in the transform matrix what would result in lossless KLT (!). Although this is not the case, we are still able to evaluate the distortion and flexibly control the coder in order to achieve least errors in the image reconstruction.

Adaptive Wavelet

3

perspective information and communication technologies” of MŠMT of the Czech Republic. REFERENCES [1] Bernas M., Páta P., Hudec R., Rezek T., Lossless and Lossy compression of Images from the OMC Experiment of Integral Project, Astrophysical Letters and Communications, pp. 429 - 432, Gordon and Breach Science Publishers, Amsterdam, 2000. [2] Burst Observer and Optical Transient Exploring System, http://www.laeff.esa.es/BOOTES/ing/index.html, 2005. [3] Castro-Tirado, A. J., et al., BOOTES: el primer observatorio robótico de España Jornadas Científicas: 250 Años de Astronomía en España Real Observatorio de la Armada, San Fernando (Cádiz), 18 Feb 2004, in press. [4] C-F-H telescope Corporation, http://www.cfht.hawaii.edu/, 2005. [5] de Ugarte Postigo, A., Mateo Sanguino, T. J., Castro Cerón, J. M., Páta, P., Bernas, M., et al., Recent Developments in the BOOTES Experimetn, In AIP Conference Proceedings 662. Cambridge: Massachusetts Institute of Technology, 553-555, 2003. [6] Effros, M., Feng, F., Zeger, K., Suboptimality of the Karhunen-Loève transform for transform coding, IEEE Transactions on Information Theory, vol. 50, Aug. 2004. [7] Hanzlik, P., Páta, P., Behind the structure of video VQ-coder, Radioengineering, vol. 13, pp. 26-32, 2004. [8] Huang, J., Mumford, D., Statistics of natural images and models, Proceedings of IEEE Conference Computer Vision and Pattern Recognition, I:541–547, 1999 [9] ISO/IEC 15444-1:2000, Jpeg2000 Image Coding System, [Konline], Dec. 2000, http://www.jpeg.org/FCD15444-1.htm, 2005. [10] Mallat, S., Wavelet Tour of Signal Processing, Academic press, 2001. [11] Páta, P., Bernas, M., In: Gamma Ray Bursts, Properties of KarhunenLoeve Expansion of Astronomical Images in Comparison with Other Integral Transforms, AIP Conference Proceedings Woodbury: American Institute of Physics, 882-886, 2000. [12] Schindler, J., Páta, P, Spatially Adaptive DWT for Image Compression, PHOTONICS PRAGUE 2005, SPIE Proceeding, 2005, in press. [13] Starck, J. L., Murtagh, F., Astronomical Image and Data Analysis, Springer, 2002.

A part of this research work has been partially supported by the research program No MSM 6840770014 “Research of

351