A fast point-pattern matching algorithm based on statistical method

TL;DR

This paper introduces a rapid pattern-matching algorithm for CCD images and stellar catalogues, utilizing a statistical method to significantly reduce computational complexity and enable millisecond-level matching.

Contribution

The paper presents a novel statistical approach for star pattern matching that outperforms traditional triangle methods in speed and efficiency.

Findings

Matching process takes only several milliseconds.

Method reduces computational complexity compared to triangle method.

Effective for images from Yunnan observatory 1-m telescope.

Abstract

We propose a new pattern-matching algorithm for matching CCD images to a stellar catalogue based statistical method in this paper. The method of constructing star pairs can greatly reduce the computational complexity compared with the triangle method. We use a subsample of the brightest objects from the image and reference catalogue, and then find a coordinate transformation between the image and reference catalogue based on the statistical information of star pairs. Then all the objects are matched based on the initial plate solution. The matching process can be accomplished in several milliseconds for the observed images taken by Yunnan observatory 1-m telescope.

Tables

Table 1. Table 1: Specifications of 1-m Telescope and CCD Detector.

Approximate focal length	$1330cm$
F-ration	$13$
Diameter of primary mirror	$100cm$
CCD field of view	7.1 ${}_{}{}^{{}^{\prime }}\times$ 7.1${}^{{}^{\prime }}$
Size of CCD array	$2048\times 2048$
Size of pixel	$13.5\mu \times 13.5\mu$
Approximate angular extent per pixel	${0.21}^{{}^{\prime \prime }}$

Equations

$$\{\begin{array}{rlr}\xi & =& \frac{\cos \delta \sin (\alpha -{\alpha }_0)}{\sin \delta \sin {\delta }_0+\cos \delta \cos {\delta }_0\cos (\alpha -{\alpha }_0)}\\ \eta & =& \frac{\sin \delta \cos {\delta }_0-\cos \delta \sin {\delta }_0\cos (\alpha -{\alpha }_0)}{\sin \delta \sin {\delta }_0+\cos \delta \cos {\delta }_0\cos (\alpha -{\alpha }_0)}\end{array}$$

$$\rho =\frac{\wp }{F}$$

$$I_{ij}=\rho \times \sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}$$

$$R_{pq}=\sqrt{(x_p-x_q{)}^2+(y_p-y_q{)}^2}$$

$$\mid I_{ij}-R_{pq}\mid <\delta F/F$$

$$\delta \theta ={\theta }_I-{\theta }_R$$

$$\mid \delta \theta -\delta {\theta }_1\mid <1^{\circ }$$

$$\{\begin{array}{rl}\xi & =\rho \cos \varphi x-\rho \sin \varphi x+c\\ \eta & =\rho \sin \varphi x+\rho \cos \varphi y+d\end{array}$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

A fast and robust pattern-matching method is very important for the modern astronomical observation and spacecraft attitude control based on star sensor. Especially for the course of the wide-field CCD imaging survey or high speed image acquisition detector such as EMCCD and SCMOS camera applied in the astronomical observations. Matching images to a stellar catalogue is a necessary prerequisite for determining the astrometric plate solution. Groth [1986] had introduced the technique of looking for similar triangles in two catalogues which used the property of affine invariant under translation, rotation, magnification and inversion between two triangles. Valdes et al. [1995] used the ratio of sides as the invariant under the coordinate transformation to define its location in triangle space and accelerated the search process by pre-sorting the triangles. Pál & Bakos [2006] described another method by culling the optimised triangles to accelerate the process. Murtagh [1992] proposed a method by characterizing a set of coordinates couples to decrease the computational complexity to $\mathcal{O}(N_{1}^{2}N_{2}^{2})$. [Mortari & Neta, 2000] proposed a k-vector search method, a fast search algorithm based on monotone function and independent of data size, to reduce the matching time.Tabur [2007] proposed two fast matching algorithms based on triangle space and vector feature. The matching speed was accelerated by matching the sparse triangles with longest side and the shortest side, and the method based on vector feature reduced the computational complexity by constructing star pairs. Heyl [2013] applied the $k-d$ method to solve the more complicated quadrilateral search problem and decreased the computational complexity dramatically from $\mathcal{O}(N_{1}^{4}N_{2}^{4})$ to $\mathcal{O}[(N_{1}^{4}+N_{2}^{4})logN_{2}]$. The computational complexity of our matching method based on statistical approach should be ${O}[(N_{1}^{2}+N_{2}^{2})logN_{2}]$ as the $k-d$ method and reduced the amount of calculation by culling the pre-sorting stars pairs. In practice, a priori knowledge of the telescope’s focal length and the CCD detector’s pixel size is common space, so the scale of the image is approximately known. Like the optimistic pattern matching described by Tabur [2007], our high-speed matching process is accomplished under a slight loss of generality.

In section 2, we give an overview of the matching method based on the statistical result of star pairs. In section 3, we analyze the performance of the method. We give our conclusions in section 4.

2 Method

2.1 Construction of Star Pairs

The first step of matching a image to a catalogue is to construct a list of sources ordered by gray values from the stellar detection and centring routine. The flat and bias corrections of the images are done before stellar detection and centring. We don’t take all of the stars in the list to construct the star pairs.We only select the brightest $m$ stars from the image source list, denoted as $\mathscr{I}$. Similarly, we sort the stars from the catalogue by magnitude and extract the $n$ brightest stars, denoted as $\mathscr{R}$. For m stars, $m(m-1)/2$ different star pairs could be obtained by means of permutations and combinations. The equatorial coordinates of catalogue’s stars should be converted to the tangent plane coordinates by the equation of 1 [Green & Green, 1985].

[TABLE]

where $\alpha$ and $\delta$ are the equatorial coordinates of catalogue’s stars, $\alpha_{0}$ and $\delta_{0}$ are the coordinates of tangent plane’s center point, $\xi$ and $\eta$ are the ideal coordinates in the tangent plane.

2.2 Searching for similar stars of the two lists

The focal length of the telescope denoted as $F$, and the pixel size of CCD detector denoted as $\wp$. The scaling relations between the ideal coordinates of tangent plane and the measured coordinate of the image can be shown as:

[TABLE]

If the telescope adopt the binning mode, the pixel size should be changed accordingly. The distance of star pair from the $\mathscr{I}$ list can be expressed as:

[TABLE]

And the distance of star pair from the $\mathscr{R}$ list can be expressed as:

[TABLE]

Then find the star pairs from the two lists that satisfy the following criteria:

[TABLE]

where $\delta F$ is the uncertainty of the telescope’s focal length. In practice we should increase the threshold appropriately. Because it’s hard to obtain the accurate uncertainty of the focal length in real time, and the scaling relations $\rho$ contains centring errors, the tangent plane’s center point error and so on. Each $I_{ij}$ might correspond to many $R_{ij}$.

2.3 Find the rotation angle between the image and catalogue

For each eligible star pair, we calculate the angle $\theta$, defined in the figure 1 :

In figure 1, the brighter star is marked as A and the other one is marked as B. The Angle is measured relative to north(although this is arbitrary) from $0^{\circ}$ to $360^{\circ}$.

The rotation angle between the coordinates of tangent plane in celestial sphere and coordinates of image the is the difference of $\theta_{I}$ and $\theta_{R}$

[TABLE]

If the $\delta\theta$ less than zero, we add $360^{\circ}$ to make the $\delta\theta$ in $0^{\circ}\sim 360^{\circ}$. We divide all the $\delta\theta$ into 360 bins($0^{\circ}\sim 360^{\circ}$), and count the amount of each bin. Looking for the most frequent of angle bins equals to find the mode of statistical result of $\delta\theta$, this value is denoted as an initial value of rotation angle $\delta\theta_{1}$. This angle is a estimated value. In order to obtain more precise angel, we need to find the subset of rotation angles $\delta\theta$ that are less then $1^{\circ}$ away from the initial angle $\delta\theta_{1}$.

[TABLE]

The angels larger than 3 times standard deviations from the mean were excluded. Notice that the final mean angle of subset, denoted as $\delta\theta_{2}$, is more accurate than $\delta\theta_{1}$. but $\delta\theta_{2}$ is not the most precise rotation angle of the affine transformation.

2.4 Determine the rough astrometric plate solution

The affine transformation between the measured coordinates of image and ideal coordinates of tangent plane in celestial sphere could be roughly expressed by a four-constant model as follows:

[TABLE]

where $\rho$ is the scale relations between two coordinate systems, and $\phi$ is the rotation angle. In this case, the $\delta\theta_{2}$ is the adopted value of the rotation angle $\phi$. Applying the equation of 3, we can obtain the translation of c and d. They are the average values after removing the outliers. Up to now, we get all the parameters of affine transformation.

2.5 Marching for all of the stars and final verification

Based on the initial plate solution obtained above, all stars detected from the image are transformed to equatorial coordinates using Equation 1 and compare to the stars of reference catalogue to find their closest match. All possible star matching pairs are reduced by the least square under a four-constant or a six-constant model. A tolerance of 2.5$\sigma$ is used, the $\sigma$ is the residual of the least square. We then use the new plate solution to do the operation above again. The repeated operation are used to avoid some false matches caused by the imprecise initial plate solution in the last step of subsection 2.4. We can use the residual of the least square to judge whether the image is successfully matched. Stars larger than 3 residual are considered to be mismatches.

3 Experimental

Several images are used to verify the reliability and speed of the method. The images are taken by the 1m telescope of Yunnan observatory. More details about the telescope and the CCD detector are given in Table 1. All of the images are taken by the $2*2$ binning acquisition mode.

The images are divided into two categories. One type of images contains asteroid and another contains M23 cluster. Reference stars are selected from the section of catalogue brighter than 18 magnitude of Gaia DR2[Lindegren et al., 2018]. Figure 2 is a successfully matched image. The green circles in the image is identified successfully, and the red ones are failed. The size of the circles is related to the magnitude of stars. The total number of stars in the catalogue is 217. Totally 77 stars are detected from the image, and 75 ones are matched to the catalogue stars successfully. Two of the failed stars in the image include an asteroid and a close binary.

3.1 Robustness of the method

Figure 3 and Figure 4 show the statistical results of rotation angels based on different number of stars from the image and catalogue. Obviously, the number of stars selected from the image m and the number stars selected from the catalogue n affects the speed of matching directly. Will it affects the success rate of matching?

We select 10 stars from the the image and 30 from the catalogue to analyze the statistical result of rotation angle $\delta\theta$ derived from the subsection of 2.3. From the Figures of 3 and 4, it can be seen that there exists a significant peak near $180^{\circ}$. The increase of both m and n like Figures 3 and 4, or the decrease of m and n like Figure 3 and 4 would not affect the correct matching. There is still a significant spike around $180^{\circ}$. In the extreme case, such as that m is 20 times n, it still match correctly even though the area of the catalog is only about twice as large as the CCD field of view.

A difficult and exhaustive matching process is to search the right matching star in a dense star field. Figure 5 shows the positions of the stars in the image of cluster M23. The image contains 135 stars, and the Gaia DR2 catalogue brighter than 18 magnitude around the center of M23 contains 33520 stars.

3.2 Time consumption analysis

The whole matching process can be divided into 3 steps. The first step is to construct the star pairs for calculating the angle distance and rotation angle of each star pair. The second step is to search the eligible star pairs from the two lists for determining the rough astrometric plate solution. The time required for two steps is independent to the total number of stars in the image and catalogue, and only relevant to the list sizes of $\mathscr{I}$ and $\mathscr{R}$. Separate timers are used to measure the performance of Step1 and Step2. Figure 6 and Figure 6 plot the time costs of Step1 and Step2 in the cases of $m=5$ and $m=10$ respectively. The time costs are the average results over 1000 trials.

The combination of $m=5$ and $n=30$ can match the most images correctly. For the non-optimal condition of partly-overlapping caused by the pointing errors or some non-stellar targets, the combination of $m=10$ and $n=50$ is used. and this combination can provide both high reliability and good performance[Tabur, 2007]. The first two steps can be completed in 3ms.

The third step is to derive the exact plate solution and match all the stars from the image to catalogue. If there are M stars in the image and N star in the catalogue, the process of finding the nearest neighbour from the catalogue of step 3 can be sped up from ${O}(MN)$ to ${O}[(M+N)logM]$[Heyl, 2013, Bentley, 1975].

4 Conclusion

This paper has outlined a efficient algorithm for matching images to a catalogue. For the most images, the star matching process can be completed in several millisecond. This method works at a slight loss of generality. We need to know the approximate focal length in advance, and attempt to match the field iteratively. It could be a viable solution when the focal length unknowns. Fortunately, the focal length of telescope and the pixel size of CCD detector are known in most cases, and the change is very little. We can magnify the uncertainty of focal length to mend this matters. This method can deal with the geometric properties about translation, rotation, local distortion and extra/missing stars. But it should be noted that this method can not handle the case of image flipping, we must know the case and take countermeasures in advance.

Ackonwledgements

We acknowledge the support of the staff of the 1-m telescope at Yunnan Observatory. This work is financially supported by the National Nature Science Foundation of China(grant nos.11503083,nos.11403101).This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.