Point Spread Functions in Identification of Astronomical Objects from Poisson Noised Image

This article deals with modeling of astronomical objects, which is one of the most fundamental topics in astronomical science. Introduction part is focused on problem description and used methods. Point Spread Function Modeling part deals with description of basic models used in astronomical photometry and further on introduction of more sophisticated models such as combinations of interference, turbulence, focusing, etc. This paper also contains a way of objective function definition based on the knowledge of Poisson distributed noise, which is included in astronomical data. The proposed methods are further applied to real astronomical data.


Introduction
Astronomy is a natural science that deals with the study of objects such as planets, moons, stars, etc. Modeling of the mentioned objects is one of the most fundamental topics in astronomical science and the branch dealing with object identification is called astronomical photometry.Under the term modeling can we understand finding of exact information about objects parameters such like its coordinates, magnitude, etc.
Previously mentioned objects are far away from the Earth and thus they appear as a bright point on the night sky.When an astronomical image, Fig. 1, is acquired, the situation is quite different and the bright point has changed and as a result the image will appear as smeared pattern.This is caused by passing of original information through the imaging system [1,2] used for an image acquisition.The result shape of captured objects is given by the impulse response of this system, also called Point Spread Function (PSF).PSF of applied imaging system is influenced by many factors and is composed as a convolution of particular PSFs of system's components.The knowledge or a good estimate of the resulting system's PSF plays a key role in astronomical photometry [3,4], when we want to know accurate information about the observed objects.
where f(k, l) are the data and n(k, l) represents noise called the dark current.This type of noise is caused by thermally generated charge, due to the long exposure times.Dark current can be simply removed by a dark frame, which maps mentioned thermally generated charge in CCD sensor.It can be considered that this type of noise is Poisson distributed [5,6,7] in the following way This article is intended to present commonly used PSF models such as Gauss and Moffat which are popular in PSF modelling are easy to evaluate and provide good approximation of analyzed objects [7,8,9,10].Another physical phenomena known as interference, focusing and atmospheric turbulence are rarely used [8,11,12,13].The main aim of this article is to compare commonly used methods with PSF's based on their combinations to show their advantage in PSF modelling.
Optimization of described PSF models is based on authors previous work dealing with detection [5] and identification [14] of analyzed objects based on analysis of dark and light frames using hypotheses testing.One of the key roles plays also the way of the objective function estimation.Objective functions determination based on hypothesis of Poisson noise is described in Section 3 as well as brief introduction of Harmony Search algorithm which was used for purpose of objective function optimization.

Point Spread Function Modeling
Astronomical photometry based on the twodimensional fitting uses the hypothesis that the profiles of astronomical point sources which are imaged on twodimensional arrays are commonly referred to as PSF [2,15], where * is a convolution operator and x, object, PSF are 2D functions, which represent the result image, the original object and the system response, respectively.PSFs can be modeled by a number of simple or more complex mathematical functions that are derived from deeper knowledge of studied problem.

Basic PSF Models
Statistical models based on different PSFs are commonly used for objects localization in astronomical science.There are usually applied two simple models, i.e., the first one is two-dimensional Gaussian function [3] where A is amplitude, x 0 , y 0 are shifts in the x − y plane, σ > 0 is its standard deviation and k, l are pixel indices as coordinates.
The second model is statistical model described by Moffat [3,16], which is a generalization of Cauchy distribution where β is a shape parameter of PDF satisfying 0 ≤ β ≤ 50.
As mentioned above, these two presented models are commonly used in astronomical photometry using (3), but do not comprise some important facts.If we consider that the light passes through the optical system before incidence onto the image sensor, than it is possible to make an approximation of the result system PSF by diffraction of circular aperture [17] where I 0 is the maximum intensity of the pattern, J 1 is Bessel function of the first order, k = 2π/λ, λ is the wavenumber, a is the radius of the aperture and θ is the angle of observation.
From physics, it is known that the diffraction phenomenon described by ( 6) is accompanied by the interference phenomenon [17] and thus we can call this relation as interference model, which is in the frequency domain expressed as where Ω > 0 is frequency and therefore adequate PSF is where H is Hankel transform [18].
Another important factor that influences resulting image is passing of the information through the atmospheric conditions.Thus, important phenomenon that influences result image is atmospheric turbulence.The interference phenomenon occurs in any optical system with constrained view.The PSF of interference is the well known Airy disc [19].The main question is whether the interference dominates the other processes, only affects them or is dominated by them.According to [19] we can express frequency spectrum of the turbulence as S(ω) = σ 2 2ψ π where ψ is the scale of the turbulence and σ > 0 represents the standard deviation of the velocity disturbance.
The role of atmospheric turbulence increases mainly in the case of wrong weather and observation conditions.
The last phenomena described in this article and that can influence result PSF is focusing [20,21], which can be expressed as follows where ρ > 0 is focusing radius and The focusing can be easily eliminated by setting ρ → 0+ in the case of perfect focusing of real telescope.But the importance of non-perfect focusing is in its combination with interference in real telescope as will be explained in the next section.

Combined Models
In this article, there are compared previously mentioned models, i.e., Gauss, Moffat, Interference, Turbulence and Focusing with more sophisticated models given by their convolutions, excluding Moffat.Thus, we can find 7 possible combinations we can use.Possible combinations are listed below and always mentioned in frequency domain with respect to the convolution theorem, as multiplying of Fourier images.
Those are Interference in combination with Turbulence (IT) which is model of constrained observation in the case of wrong atmospheric condition.Fortunately, the effect of turbulence can eliminate Airy disc around every point light source.
Interference plus Focusing (IF) which is model of constrained observation in the case of wrong focusing.But this non-perfect focusing is very pragmatic way of Airy pattern suppression.Therefore, the IF model is the favorite candidate to modeling of images from real but unknown telescopes.

Gauss model in combination with Interference (GI)
where F G is Fourier transform of (4) and can be expressed as This model is similar to IF but less pragmatic and used only for completeness.
Next models are the other useful combinations of basic ones.Those are Gauss with turbulence (GT) as Focusing with Turbulence (FT) Gauss, Interference and Focusing (GIF) and finally combination of Interference, Focusing and Turbulence (IFT) This approach was used due to the reason that analytical convolution of S(ω) to s(r) is impossible, thus the data analyzed by combined models are processed only in the frequency domain, subsequently transformed into spatial domain and properly modified for amplitude A, x 0 and y 0 shifts optimization.
Parameters of particular single models and their lower and upper bounds can be found it Tabs. 1   Tab. 2. Lower and upper bounds of single models parameters.

Objective Function Definition
Optimization of models introduced in Sec.2.1 is based on minimization of objective function.Its inference can be performed using Least Square Method (LSM), but as mentioned in [5], data acquired by astronomical CCD camera are Poisson distributed.Thus it is upon the place to use different approach based on Maximum Likelihood Estimate (MLE), see [14].
In statistics, MLE is a method of estimating the parameters of a statistical model (4).When it is applied to a data set and given statistical model, MLE provides estimates for the model's parameters.For a fixed set of data and certain statistical model, it produces a distribution that gives to the measured data the greatest probability, i.e., estimated parameters maximize the likelihood function [22,23].
Let us now consider astronomical image defined by (1), noise model (2) and average dark frame Model image with astronomical objects described by PSF model can be derived from (1) by replacing expression f(k, l) in the following way where (k, l) ∈ D M×N , M and N are dimensions of the rectangle region of interest D and f(k, l, p) is PSF model of astronomical object with vector of parameters p.
When it is supposed that the data x are Poisson distributed with number of occurrences λ then it is possible to write that When the MLE is used to (24) and x is replaced by (22), then the opposite likelihood function can be written as where x(k, l) is the analyzed light image, d(k, l) presents appropriate average dark frame and f(k, l, p) is the diffusion model, whereof parameters are estimated.
Combination of ( 24) and ( 25) leads to the final form of function φ where c is some constant.The constant c is only data depending and can be set to satisfy φ ≥ 0 and obtain For the purpose of PSF modeling, astronomical objects were classified into three classes based on the bit depth of analyzed image.Processed data were acquired in the 16 bit depth, thus the maximum intensity is 65,535.Intervals of intensity values were uniformly divided into three classes, which can be written as follows: • small object-maximum intensity in the analyzed area is less than 21,845, • medium object-maximum intensity in the analyzed area is higher than 21,845 and less than 43,690, • large object-maximum intensity in the analyzed area exceeds 43,690 and the top is given by the system resolution properties, thus 65,535.
Chosen objects that were used for an application of proposed methods can be seen in Fig. 2.

Objective Function Optimization
For the purpose of objective function optimization, Harmony Search algorithm was used due to its good results related to estimated quality and variability of results based on authors' previous research [24].
Harmony Search (HS) [25,26] is inspired by musician improvisation process.It imitates the natural phenomenon of musicians' behavior when they cooperate the pitches of their instruments together to achieve a fantastic harmony as measured by aesthetic standards.This musicians' prolonged and intense process led them to the perfect state.It is a very successful metaheuristic algorithm, which can explore the search space of a given data in parallel optimization environment, where each solution (harmony) vector is generated by intelligently exploring and exploiting a search space.

Results
Tabs. 3-8 summarize results of single and combined models optimization, where there are mentioned minimum, maximum, average values of φ values and standard deviation.There are also mentioned parameters of best object model estimate gained by HS algorithm used in this article.
Figs. 3 and 5 shows graphical representation of φ min for small, medium, and large objects.The best three basic mod-els are highlighted in Tabs. 3 -5.As can be observed, traditional Gaussian and Moffat PSF models were over-performed by interference model in all cases.But in the case of large object, the focusing model is better than the interference one.As seen in the case of basic models, the phenomenon of interference in the telescope dominated.The best three combined models are highlighted in Tabs.6 -8.The bests of them were better in PSF fitting than the best individual filters in all cases, which is the main argument for the mixed model application.The best models were GI, IF, and GIF but in various order which depended on the object size.Therefore, the combination of the interference with Gaussian model and/or focusing offered the best results.
The other characteristics (φ max , φ avg , φ std ) included in Tabs. 3 -8 are rather useful for the evaluation of optimization process complexity and reliability than for optimum model selection.From φ std point of view, the models consisting of Gaussian, Moffat or Turbulence sub-model are easy to fit via several runs of Harmony search.The evaluation of φ max , φ avg comes to the same result.

Conclusion
From the real perspective, it is obvious that application of combined models is better than use of single models.They include multiple influences on real information passing through atmosphere and optical system of used astronomical equipment.
We found that IF, GI and GIF are the most suitable models.If we look on models, where there was used turbulence of atmosphere, those were always worse than the three aforementioned models.Turbulence itself reached worse results in case of simple models application thus from modeling point of view it seems inappropriate to use it.
Therefore, any PSF model consisting interference phenomenon is recommended for the fitting of astronomical objects but repeated application of optimization heuristics is necessary for reliable fit.In contrast, Gaussian, Moffat and turbulence models are very easy to fit.