Automated Meteor Detection by All-Sky Digital Camera Systems

Abstract

We have developed a set of methods to detect meteor light traces captured by all-sky CCD cameras. Operating at small automatic observatories (stations), these cameras create a network spread over a large territory. Image data coming from these stations are merged in one central node. Since a vast amount of data is collected by the stations in a single night, robotic storage and analysis are essential to processing. The proposed methodology is adapted to data from a network of automatic stations equipped with digital fish-eye cameras and includes data capturing, preparation, pre-processing, analysis, and finally recognition of objects in time sequences. In our experiments we utilized real observed data from two stations.

1 Introduction

The astronomical origin of meteors was discovered already in the nineteenth century. Today, we look at meteors as messengers carrying information about the nearest space. However, to decipher this information a general knowledge concerning the size and frequency of meteors from individual directions is required. Therefore networks of observation stations are built covering vast territories, e.g. over Central Europe (Spurný et al. 2007). In rare cases, when a celestial body reaches the Earth surface as a meteorite, the precise computation of the meteoroid trajectory is essential to determine the impact location. Finding a meteorite on the Earth’s surface is always an extraordinary event, especially if it is located by an accurate path tracking calculation.

Meteorite falls are products of an interaction of celestial bodies (fragments of comets and asteroids) with Earth’s atmosphere. A meteorite found on Earth’s surface is always an extraordinary event, especially if it is based on precise computation of path tracking. Following chemical analysis of the meteorite found can provide another information about the composition and properties of the original (parent) bodies. A meteorite found on the basis of instrumentally documented falls is sometimes called as a meteorite associated with the lineage. The very first one is described in Ceplecha (1961).

Formerly, meteor observations were captured using traditional analogue cameras equipped with photographic material (films or plates). Exposition usually lasted all night, thus examining one snap per night visually, typically from several stations, did not pose severe problems. State-of-the-art stations with extensive camera networks are based on automatic all-sky CCD systems. Obviously, such sophisticated automated systems churning out terabytes of image data per night require a relevant data processing system.

The main aim of the proposed methodology is to automatically sort snaps containing meteor tracks at maximum confidence level. This general goal includes many partial tasks depending on observation conditions and data character sets, in particular distinguishing between false positives and false negatives in the search stage is extremely important. A false negative means a meteor was not found at all or the cause of the light track was identified incorrectly as something else. Each false negative is a significant failure of the algorithm. On the other hand, a false positive means a different feature was identified as a possible meteor. Likewise, an increasing number of false positives leads to a progressive failure of the algorithm.

Fig. 1

Flowchart of the proposed procedure for one snap. The block “Tests of the lines” can be divided into five subblocks

The flowchart of the proposed methods can be seen in Fig. 1. The organization of the paper follows the flowchart. Section 3 deals with data acquisition, Sect. 4 describes pre-processing (Masking of ground objects, difference adjacent frames, and hot pixel removal), Sect. 5 explains the Hough transformation with and without the fish-eye modification, and Sect. 6 discusses tests of the found lines. Section 7 shows our results, Sect. 8 is a brief note about computing complexity and Sect. 9 concludes the paper.

1.1 State of the Art

There are several observation networks for meteor search across the world. The IMO Video Meteor Network was originally established in Germany. By December 2013, the network has grown to 88 cameras operated by 49 observers in 16 countries, see Molau (2001). MetRec is based on video cameras with the standard frame rate 25 or 30 fps, then a meteor trace is captured in a number of images combined in a video sequence.

The continuosly growing SonotaCo network, which consisted of 70 observation stations in the year 2010, operates in Japan. The system starts recording only once a moving object appears and utilizes the UFOCapture, UFOAnalyzer, and UFOOrbit software. General information on the network is available on its web pages (SonotaCo 2011). It uses the star catalog for the field-of-view alignment and the meteor stream catalog to specify meteor detection. We believe in advantages of subtraction of adjacent frames for star removing.

The Californian CAMS (Jenniskens et al. 2011) consists of three observatories each operating 60 video cameras. They monitor the sky above $31^\circ$ elevation. The network utilizes its own software which includes modules from other packages, e.g. from MeteorScan. The trail search is based on the Hough transformation.

The Polish Fireball Network (PFN) consisted in the year 2016 of 36 continuously active stations with 57 sensitive analogue video cameras and 7 high resolution digital cameras (PFN 2004; Wiśniewski et al. 2017). A special software package “PyFN” for trajectory and orbit determination was developed for this net. The meteor trails are searched by MetRec (Żołądek et al. 2007).

The European viDeo MeteOr Network Database (EDMOND) is a database of video meteor orbits resulting from cooperation and data sharing among several European national networks and the International Meteor Organization Video Meteor Network, IMO VMN, see Kornoš et al. (2013). The corresponding network is called EDMONd. At present it consists of observers from the following national networks: BOAM (France); BosNet (Bosnia); CEMeNt (cross-border network of Czech and Slovak amateur observers); CMN (Croatia); FMA (Switzerland); HMN (Hungary); IMO VMN; MeteorsUA (Ukraine); IMTN (Italy); NEMETODE (United Kingdom); PFN (Poland); Stjerneskud (Denmark); SVMN (Slovakia); and UKMON (United Kingdom). It utilizes the UFOOrbit software from SonotaCo while collecting data from 155 stations.

The SPanish Meteor Network (SPMN) consists of similar low-scan-rate all-sky CCD cameras as the AFO4 (Autonomous Fireball Observatories) Czechia network (Spurný et al. 2007). The software (Trigo-Rodríguez et al. 2008) computes meteor trajectories and orbital data. It uses a modification to the traditional Hough transformation called “a phase-coded disk”. The ground objects are removed by subtraction of adjacent frames; the slowly moving stars create a typical pattern of one positive and one negative crescents removed by a special procedure.

References to other networks can be found in Wiśniewski et al. (2017); the good summary is also in Gural (2008). The software cited in part inspired our work, however we could not use it directly as it was developed for a different type of capturing devices, typically for video cameras. In some cases, we decided to use an alternative approach to the existing methods, e.g. to subtract images with correction of Earth’s rotation instead of non-corrected images.

Likewise, one must mention a diploma thesis (Fajfr 2013) that processed data from the first digital camera placed at the Ondřejov station. The student implemented the idea of searching for linear objects in an image using the Canny edge detector.

2 Source of Data

We tested our operating procedures on the real data from the AFO network built and upgraded to Digital Autonomous Fireball Observatories (DAFO) thanks to Pavel Spurný and his team in the Czech Republic; details can be found in Spurný et al. (2007). DAFO consists of 11 stations distributed over Czechia distanced less than 200 km from each other. In addition to the necessary electronics each autonomous station consists of two special fish-eye digital cameras observing the entire sky and operating in a dual mode.

2.1 Digital Camera Parameters

Each camera alternates 35 s of capturing and a 25 s pause with an overlap of 5 s. The active phase of the camera is interrupted with a frequency of 15 Hz and a duty cycle 50%, i.e. 1/30 s capturing and 1/30 s pause; each 15th pause is omitted. This regular sampling is an electronic replacement for the previously used “rotating shutter” providing time-marks in the meteor trail. It is based on an LCD component alternating transparency and opacity. Using these parameters we can estimate the velocity of a flying object.

Frequencies of the observation mode are set so that the light meteor trail leaves a regularly interrupted curve, while the slower moving objects evince the trail either interrupted irregularly as airplanes or without interruptions as satellites. The size of the color 8-bit snaps is $5472\times 3648$ pixels. Color significance is minor as meteor’s signs are white, therefore we work with brightness only

$$\begin{aligned} Y = 0.299 R + 0.587 G + 0.114 B. \end{aligned}$$

(1)

A captured object that is too colorful typically signifies the trail originator not being a meteor. Such colorful light trails are rare, therefore we do not test colorfulness in the image. Available image data for our experiments are from two stations, Ondřejov and Churáňov. An example of the Ondřejov’s original is in Fig. 2.

Fig. 2

Example of the Ondřejov station data, snap from 12th January 2014, 17:28:30, camera 0

2.2 Fish-Eye Camera Model

Fish-eye cameras are used in AFO; their advantage is clear: one camera can observe the whole sky. The drawback is that the originally straight trails of flying objects are displayed as arcs requiring error compensation. We introduce the fish-eye camera simple model as projection on spherical surface, see Fig. 3.

Fig. 3

Simplified model of the fish-eye camera

A pixel distant $\varrho$ from the optical axis is converted to a pixel distant $\varrho '$ from the optical axis of a traditional pin-hole camera. We can easily derive $\varrho = R \alpha$ ($\alpha$ in radians) and $\varrho ' = R \tan \alpha$, from that

$$\begin{aligned} \varrho ' = R \tan \left( \frac{\varrho }{R}\right) \end{aligned}$$

(2)

for conversion from fish-eye to pin-hole and

$$\begin{aligned} \varrho = R \arctan \left( \frac{\varrho '}{R}\right) \end{aligned}$$

(3)

for conversion from a pin-hole to a fish-eye cameras. The sphere radius R was experimentally estimated as 0.3 of the smaller image size 3648, i.e. 1094.4 pixels. Due to model imprecision we estimate an error between 2 and 3 pixels, which cannot be removed by an improved estimate of R.

The entire conversion procedure begins with conversion to polar coordinates

$$\begin{aligned} \vartheta& {}= \displaystyle \arctan \left( \frac{y-y_c}{x-x_c}\right) \\ \varrho& {}= \sqrt{(x-x_c)^2+(y-y_c)^2}, \end{aligned}$$

(4)

then $\varrho$ is converted to $\varrho '$ by (2) and back to Cartesian coordinates

$$\begin{aligned} x'& {}= x_c + \varrho '\cos \vartheta \\ y'& {}= y_c + \varrho '\sin \vartheta . \end{aligned}$$

(5)

The back conversion is analogous (3) assuming the optical center {$x_c$,$y_c$} is at the real center of the image, i.e. $x_c=2735.5$, $y_c=1823.5$ in zero-based coordinates.

The pin-hole camera image is significantly larger compared to the one obtained using a fish-eye camera, theoretically $\varrho =\pi /2$ is mapped to infinity, when the whole hemisphere is used. These extreme positions, however, are not used in practice because the real horizon is higher. If possible, computation of the whole pin-hole camera image should be avoided due to extreme memory demands and higher computing complexity of certain procedures, as illustrated further.

3 Bright Tracks on the Snaps

Images collected from the AFO network incorporate various other light tracks not caused by meteors. These tracks originate from both terrestrial objects at the circle of the horizon and objects in the inner part of the image, e.g. the Moon, stars, planets, airplanes, satellites, and clouds. Recognition of objects other than meteors and their exclusion from the search process is the most demanding phase of the proposed methodology principally as the same class of objects (e.g. class of airplanes) may produce differing time variable bright trails. Examples of the light tracks of satellites and airplanes are in Fig. 4.

Fig. 4

Light tracks of satellites and airplanes: a ISS, b other satellite, c airplane, d other airplane

3.1 Ground Objects

The ground objects in the image are sited near the horizon and are typically stable overnight. In principle, these objects can be divided into illuminated (any reflective surfaces, windows of buildings, etc.) and non-illuminated (trees, poles).

We propose a binary image as a mask to remove these objects which is a very reliable approach, however, with one drawback. The mask for each camera in the pair ought to be created separately with corrections introduced whenever the surrounding changes. The non-illuminated objects are darker than the sky and can be segmented by region growing. The mask yields improved results when the background is dilated by several pixels, i.e. the remaining sky is eroded. We use the morphological operation with circular structure element 5 pixels in diameter resulting in filling the holes both in the sky and beyond. Possible remaining illuminated objects in the sky part of the mask are removed manually. An example of the mask and its application is in Fig. 5.

Fig. 5

Ground object erasing: a original snap, b mask for the Ondřejov station, camera 0, c snap without ground objects. Empty lateral margins are omitted. d–f Enlarged details of the edges corresponding to the upper frames

3.2 The Moon, Stars and Planets

Extending the observation mode to nights between the first and last quarter of the Moon introduced several problems to the image processing. The Moon, being the brightest object in the image, will hardly be recognized as a meteor, however it has other unfavorable impacts. Firstly, weaker meteors in the vicinity of the Moon are practically invisible and secondly, a cloud edge illuminated by the Moon, could be confused with a meteor.

Another problem we have to take into account are stars, whose positions create a straight line in a certain section. At first glance, such features indeed appear to look like a meteor. These collinear group of stars are present in the sky throughout the whole night. Search algorithms then filter out every snap of the night as “containing meteor”, i.e. false positive. Although minor irregularities can be seen in the distances of the adjacent stars, we were not able to propose a sufficiently reliable frequency test that would exclude the star lines from the list of meteors.

The solution is based on the removal of all stable objects from the snaps. In a snap sequence of the whole night we subtract the adjacent snaps from the current one, hence removing stable objects and preserving the ones moving. However, this process induces a complication associated with the Earth’s rotation. A 1-min difference between adjacent snaps from the same camera registers a quarter of a degree Earth rotation. This rotational angle represents up to 12 pixels in our image data and cannot be neglected. To compensate for the rotation would be trivial in snaps from a pin-hole camera, but for the fish-eye is more complicated.

The coordinates of each pixel are converted to the pin-hole by the procedure described in Sect. 2.2, then are rotated around the celestial north pole {$x_o$,$y_o$} by the angle $\beta$ and finally inverted to the fish-eye coordinates. The rotation around the north pole is performed using the formulas

$$\begin{aligned} x'& {}= x_o + (x-x_o)\cos \beta - (y-y_o)\sin \beta \\ y'& {}= y_o + (x-x_o)\sin \beta + (y-y_o)\cos \beta , \end{aligned}$$

(6)

where $\beta =0.25^{\circ }$ for the previous snap and $\beta =-\,0.25^{\circ }$ for the next snap.

A separate task is the location of the celestial north pole {$x_o$,$y_o$} in the snaps. As it is near Polaris, we defined a square, where the only significant star present is Polaris. The maximum brightness pixel detected within this square represents the center of Polaris, see Fig. 6. An elliptical regression (Fitzgibbon et al. 1996) of the Polaris trail is then computed with the ellipse center representing the coordinates of the celestial north pole. This procedure is repeated for each camera separately, see Table 1.

Fig. 6

Polaris positions during the observing night including the fitted ellipse

Table 1 Zero-based celestial north pole coordinates

Since our model of the fish-eye camera is not perfect, the stars would not be subtracted completely. This is why we use a “maximum” type filter for both previous and next snaps. This filter substitutes the value of each pixel by the maximum from its neighborhood in the form of a circular disk 7 pixels in diameter. Then the previous and next snaps are subtracted with suitable weights from the current snap.

The overall brightness is changing during the observing night in dependence on the Moon-/Sun- -rise/-set. We need to compensate for these disturbing effects manifested throughout the entire sequence. The solution is based on subtraction regarding to contrast and brightness bias normalization. The subtracting formula is

$$\begin{aligned} f_c'(x,y)=f_c(x,y)-\left( \left( f_p(x,y)-\mu _p\right) \frac{\sigma _c}{\sigma _p} +\mu _c+\left( f_n(x,y)-\mu _n\right) \frac{\sigma _c}{\sigma _n}+\mu _c\right) /2, \end{aligned}$$

(7)

where $f_c(x,y)$ is the current snap, $f_p(x,y)$ is the previous snap and $f_n(x,y)$ is the next snap. The previous and next snaps are after Earth’s rotation compensation, ground object masking and maximum filtering. The symbols $\mu _c$, $\mu _p$ and $\mu _n$ are brightness averages of the current, previous and next snaps, respectively, $\sigma _c$, $\sigma _p$ and $\sigma _n$ are their brightness standard deviations.

At the beginning of the sequence, there is no previous snap, similarly at the end, there is no next snap. In these particular cases we subtract only one snap without the coefficient 1/2. After subtraction, the negative values are substituted by zeros, otherwise the possible light trails in the previous or next snaps could violate the search of trails in the current snap.

The movement of planets is negligible with respect to the Earth’s rotation, therefore we subtract them together with the stars and the Moon. The effect of star subtraction can be seen in Figs. 7 and 8. The Moon and the stars are off, but the clouds, meteors and other moving objects remain, see the distant airplane at the bottom of Fig. 7.

Fig. 7

Detail of meteor trail: a before and b after the subtraction of previous and next snaps. The contrast of the snap (b) has been six times increased

Fig. 8

The current snap after the subtraction of previous and next snaps. The Moon is gone, an illuminated cloud remains

4 Hot Pixels and Other Impurities

Despite the elimination of the stars performed in the previous step, several strange objects remain in the snaps. At first glance, they look like small stars. However, after examining details, we found several differences. Firstly, these stains appear on one snap only. The previous and the next snaps have stains distributed elsewhere, therefore they cannot be subtracted. Secondly, they are more symmetrical than stars which, due to the Earth’ rotation, appear as ellipses with small eccentricity. Thirdly, some objects look like small disks while others resemble the number five on a dice. The origin of these objects is unknown. The colorful ones are probably hot pixels, i.e. overexposed solitaire pixels on the CCD sensor. The white objects may be cosmic radiation or satellite reflections. Due to the unknown origin of the objects a proper terminology can hardly applied, hence we call all these impurities by the common name “hot pixels”.

In certain cases a collinear group of hot pixels can cause a false positive. These false positives are rarer than those caused by the stars, nevertheless a few such cases were presented in the tested data. So, we decided to reduce the possibility of any other false positive.

We suggested two tests how to deal with these impurities. The first test (solitariness test) supposes the distances of the hot pixels are bigger than the distances of the dashes in the meteor trail. We have two circular neighborhoods for each pixel, the inner one with a diameter of 5 pixels and the outer one with a diameter of 25 pixels. When the sum of brightness in the inner neighborhood is greater than 0.1 (we use the range of brightness from 0 to 1) and the sum of brightness in the outer neighborhood without the inner one is less than 0.4, then the inner neighborhood is erased.

The second test (eccentricity test) measures eccentricity of the characteristic ellipse of a hot pixel. The geometric moments of the image f(x, y) are defined

$$\begin{aligned} m_{pq}=\int \limits _{-\infty }^{\infty }\int \limits _{-\infty }^{\infty } (x-x_t)^p(y-y_t)^q f(x,y){\mathrm {d}}x{\mathrm {d}}y, \end{aligned}$$

(8)

where {$x_t$,$y_t$} are coordinates of some significant point in the image. Here, we sum over a circular neighborhood with a diameter of 9 pixels. {$x_t$,$y_t$} are coordinates of the neighborhood center (i.e. the current pixel). The characteristic ellipse has the same moments for $p+q$ from 0 to 2 as the original object. When the image has a local maximum and the eccentricity is less than 0.02, we erase the neighborhood. The eccentricity is

$$\begin{aligned} e_c=\frac{\sqrt{(m_{20}-m_{02})^2+4m_{11}^2}}{m_{00}^2}=\frac{a^2-b^2}{\pi ab}, \end{aligned}$$

(9)

where a, b are semi-major and semi-minor axes of the ellipse. The traditional relative eccentricity is computed as

$$\begin{aligned} e=\sqrt{\frac{2\sqrt{(m_{20}-m_{02})^2+4m_{11}^2}}{m_{20} +m_{02}+\sqrt{(m_{20}-m_{02})^2+4m_{11}^2}}}=\sqrt{1-\frac{b^2}{a^2}}, \end{aligned}$$

(10)

but we have a better experience with $e_c$. Examples of the hot pixels are in Fig. 9a, b.

Fig. 9

Hot pixels: a circular and b the number five on a dice; c, d performance of the solitariness and eccentricity tests and hot pixels erasing

5 Moving Objects and the Hough Transformation

Finally, we eliminated the static objects in the snaps (both ground and celestial) and now we address the moving objects and the pattern recognition of their trails. A familiar method for the search of linear grouping in images is the Hough transformation (HT).

Any straight line can be described by the distance r from the coordinate origin and the angle $\theta$. Traditionally, $\theta$ is the angle from the y-axis to the straight line, see Fig. 10, then its description is

$$\begin{aligned} r=x\cos \theta +y\sin \theta . \end{aligned}$$

(11)

All points {x,y} satisfying this equation for a fixed r and $\theta$, lie on the straight line.

Fig. 10

Traditional HT

The computation consists of a double outer loop over pixels {x,y} of the image and of an inner loop over $\theta$. We obtain an output array with two indices, both for r and $\theta$. r has been computed by (14) for each x, y and $\theta$. The f(x, y) supplements the output array element with indices $\theta$ and corresponding r. If the computation yields a significant local maximum in the output array, the corresponding straight line has a significant response in the image.

5.1 Fish-Eye Modification

The fish-eye camera distorts originally straight lines to arcs. One possibility would be to search them as circles with large radii, hence a modification of the HT for circle search is needed. Its computing complexity is higher, because the circle has three parameters compared to two parameters of the straight line, requiring one additional loop in the algorithm. This modification is usually used for search of circles that are fully included in the image. If our circles had their centers far outside the image, the memory and time demands would be enormous. That is why we decided to search only such arcs that would appear as straight lines in the pin-hole camera.

Our modification has the same loops over pixels and angles as the original HT. First, we convert the pixel coordinates to polar by (4), then we convert the fish-eye radius $\varrho$ to the pin-hole radius $\varrho '$ by (2). The distance of the pin-hole straight line from the coordinate origin is computed from

$$\begin{aligned} r' = \varrho '(\cos \vartheta \cos \theta +\sin \vartheta \sin \theta ). \end{aligned}$$

(12)

The pin-hole distance $r'$ is converted to the fish-eye distance r by (3), i.e.

$$\begin{aligned} r = R \arctan \left( \frac{r'}{R}\right) . \end{aligned}$$

(13)

Finally, f(x, y) is added to the output array element with indices $\theta$ and r as in the traditional version.

5.2 Parameters of the HT

The HT application is inefficient in the whole image. The step of angle increment must be very small increasing computing complexity. In spite of this, the prospect of finding short weak meteor trails is lower than in smaller parts of the image. Therefore we decided to split the image into $12\times 8$ subwindows each consisting of $456\times 456$ pixels and to compute the HT in each subwindow separately.

There is a danger that the meteor track lies on the border between the subwindows. The minute size of each individual part of the track may prevent detection in any subwindow. Therefore we use 10% overlap of the subwindows, i.e. 47 pixels. Subwindows in the center of the image then have the size of $550\times 550$ pixels. The empty subwindows in the first and last two columns are not used for computation, thus in practice we have $8\times 8$ subwindows.

After multiple experiments, we determined the optimal angle step of $\theta$ as $1.5^{\circ }$ and the shift step of r as 5 pixels. This means r is rounded to 5 pixels and ranges from $-\,r_m$ to $r_m$, where $r_m$ is the image diagonal, i.e. 2671 bins for the shift, which is much more than 120 bins for the angle. The complete size of the output array is $2671\times 120$. The 5 pixel step for r also means the 5-pixel thick line is integrated into one element of the array and the prospect of finding weak meteors is higher than for less steps.

The global maximum and significant local maxima in the output array are found. After a maximum is found, its neighborhood $59\times 9$ is cleared to guarantee sufficient distance between the maxima. Each detected maximum {r,$\theta$} corresponds to a closed curve and we must find their parts in the searched subwindow. There is a loop over an angle $\phi$, the corresponding point lying on the curve is computed

$$\begin{aligned} \rho & {}= R\arctan (R(\cos \phi \cos \theta +\sin \phi \sin \theta )/r)\nonumber \\ x & {}= x_c+\rho \cos \phi \nonumber \\ y & {}= y_c+\rho \sin \phi . \end{aligned}$$

(14)

If the point {x,y} lies in the searched subwindow, it is included in the tested trace, see thick parts of the circle in Fig. 11b.

If all the tests fail, the traditional HT for a straight lines is computed. The meteor trails passing through the image center are not distorted and can be found more reliably without the fish-eye modification. We use a different approach to search for the intersection of the found straight line with the subwindow. We have found maximum {r,$\theta$} by the HT. The zero-based coordinates of the vertices of the subwindow with size $x_e\times y_e$ are

$$\begin{aligned} x_w & {}= \{0,x_e-1,x_e-1,0,0\},\nonumber \\ y_w & {}= \{0,0,y_e-1,y_e-1,0\}. \end{aligned}$$

(15)

We parametrize the border lines of the subwindow as

$$\begin{aligned} x=x_{wi}+t_i(x_{w,i+1}-x_{wi})\\ y=y_{wi}+t_i(y_{w,i+1}-y_{wi}),\\ \end{aligned}$$

(16)

where $i=1,\ldots 4$ and $t_i$ is the parameter. If $0\le t_i\le 1$ for some i, the corresponding point lies on the border of the rectangle. The HT found the line inside the rectangle, therefore at least two out of four values $t_i$ must satisfy it (equality can hold for more than two points). The specific values can be found

$$\begin{aligned} t_i=\frac{r-x_{wi}\cos \theta -y_{wi}\sin \theta }{(x_{w,i+1}-x_{wi})\cos \theta +(y_{w,i+1}-y_{wi})\sin \theta },\quad i=1,\ldots 4. \end{aligned}$$

(17)

Let us denote the parameter values between 0 and 1 as $t_a$ and $t_b$. The final points of the found straight line segment have then the following coordinates

$$\begin{aligned} x_a & {}= x_0+x_{wa}+t_a(x_{w,a+1}-x_{wa}),\nonumber \\ y_a & {}= y_0+y_{wa}+t_a(y_{w,a+1}-y_{wa}),\nonumber \\ x_b & {}= x_0+x_{wb}+t_b(x_{w,b+1}-x_{wb}),\nonumber \\ y_b & {}= y_0+y_{wb}+t_b(y_{w,b+1}-y_{wb}), \end{aligned}$$

(18)

where {$x_0$,$y_0$} are the coordinates of the left top corner of the subwindow. The situation is drawn in Fig. 11a.

Fig. 11

Intersections of the rectangular subwindow with the found curve. a Line segment can be found as connection of two points with the curve parameter t between 0 and 1. b Analytical solution for circles is more complicated, it is easier to trace the circle pixel by pixel

Examples of meteor trails are in Fig. 12.

Fig. 12

Meteor trails found by the HT: a with the fish-eye modification, b without the modification

6 Tests of the Lines

Each line detected by the HT must be traced to find the trail followed by a test to recognize a meteor or some other object. First, the curve is reparametrized, so the distance of the adjacent points is 1 pixel (so called reparametrization by distance).

The meteor trail can begin or finish somewhere in the middle of the subwindow. We have a threshold $\theta _b=1.5\%$ for the maximum brightness. The parts of the trail near the borders under this threshold are omitted. Only the trail part between the first and the last over-threshold points is tested. The linearity test is the first.

6.1 Linearity Test

This test serves to distinguish between linear traces and planar objects (mostly clouds). A mask of pixels with a distance less than 8 pixels from any point of the curve is created. The average brightness $m_b$ in this mask is computed. Another mask is made from the boundary points of the first mask. The average brightness in this mask is $b_b$. The coefficient $m_b/b_b$ must be greater than 3, otherwise the trace is not sufficiently linear. Pixels outside the linear trace have a significantly lower brightness than those inside. Pixels in the planar object have a similar brightness everywhere. Figure 13 shows a false trace removed by this test, $m_b/b_b$ equals only 1.068.

Fig. 13

Linearity test. a Original part of the snap, where a false trace was found in clouds. b Preprocessed part of the snap—the static objects have been removed and the contrast is increased. The green curve accentuates the detected track. c The red strip delimits the curve neighborhood. Average brightness is calculated both on the border and inside the strip.

6.2 Fourier Test

Interesting experiments with a rotating shutter can be found in Bettonvil (2008), where it is demonstrated how the precision of the meteor speed estimate depends on the rotation speed of the shutter. However, we have no access to the hardware, therefore we cannot change the frequency of our “rotating shutter”.

This Fourier test should verify, if the trace contains spaces between dashes from the electronic rotating shutter. First, we look for the longest segment of the trace. The spaces in the segment with a brightness less than $\theta _b=0.015$ must not be longer than 20 pixels. Then we compute a 1D Fourier transformation of the longest segment with values $s_i$ of length d

$$\begin{aligned} a_j=\sum \limits _{i=0}^{d-1}s_i e^{-2\pi ij\sqrt{-1}/d},\quad j=0,\ldots d-1. \end{aligned}$$

(19)

The amplitude spectrum $|a_j|$ has a significant local (often global) maximum at 0 equaling the sum of brightness in the trace. If it is a meteor trace, there should be another local maximum corresponding to the rotating shutter. The spectrum is normalized

$$\begin{aligned} a_j'=100\,|a_j|\,\frac{\sqrt{j+1}}{d+1}. \end{aligned}$$

(20)

In the Fourier test, the spectrum is analyzed from 3 to $d_h=\min \{d/2,64\}$. First, we look for the first local minimum $p_n$ less than 0.9 between 3 and d. Let us suppose its position is n. Then we look for the maximum $p_m$ between n and $d_h$ in the position m. The number of local maxima between n and $d_h$ higher than 0.7 is $p_x$. The trace passes through the Fourier test, when $d>20$, $p_n/p_m<0.9$, $p_m>0.6$, $m>3$, $d/m>3$, and $1\le p_x\le 3$.

An example of the spectrum $a_j'$ used for the Fourier test is in Fig. 14. There is a significant local maximum $p_m=9.85$ at $m=32$. The other features are: $d=351$, $p_n/p_m=0.132$, $d/m=10.97$ and $p_x=1$.

Fig. 14

Normalized amplitude Fourier spectrum: a covered low frequencies only, b including higher frequencies

6.3 Duty Cycle Test

The duty cycle of the rotating shutter simulation is 50%, i.e. the opening time is one half of the whole time. The duty cycle in the trace should be similar, the percentage of dashes in the length of the trace is about 50%.

The brightness in the trace fluctuates; we need to search the point, where the dash changes to the dark space. Here, the light curve crosses a middle level between a minimum in the space and a maximum in the dash. We compute this level by filtering the light curve by a Gaussian filter of a length 21 and standard deviation 2.5 pixels. It is subtracted from the original signal and the level is shifted by 1/25 of the global maxima, see Fig. 15. The results are more stable at this shifted level.

The duty cycle c is then $s_c/d_c$, where $s_c$ is the number of pixels with the value above the level, and $d_c$ is the distance of the first and last such pixel. The signals passing through the test meet the requirements $20\% \le c \le 60\%$.

Fig. 15

Duty cycle computation: a original light curve in the trace, b subtraction of local mean values. The horizontal line shows the level, where the duty cycle is computed, $c=44.6\%$ in this case

6.4 High-Frequency Test

The spectrum of meteors should not contain a much higher frequency signal than that from the rotating shutter. On the other hand, airplane traces and clouds contain much high-frequency noise. Occurrence of the high frequencies is measured between quarter and half frequencies of the sampling (between Nyquist frequency and its half)

$$\begin{aligned} v=\frac{1}{d/2-d/4+1}\sum \limits _{j=d/4}^{d/2}a_j'. \end{aligned}$$

(21)

It is related to the value $p_m$ of the spectrum on the rotating shutter frequency m. The signal passing through this test complies $v/p_m < 0.14$. The spectrum in Fig. 14b has $v/p_m=0.032$.

7 Real Data Experiments and Results

The proposed methodology of automatic detection of meteor traces has been tested and validated using data from two stations—Ondřejov and Churáňov. In total we processed data from 8 observing nights (1 Ondřejov, 7 Churáňov), two cameras at each station. A night sequence of image data consists of approx. 1700 snaps depending on the length of observation. An overview of available data is in Table 2.

Table 2 Number of snaps

7.1 Found Meteor Trails

The thresholds of the tested data were intentionally set to avoid false negatives. All meteor trails were found in the whole sequence and none was completely missed. Only meteor trails that were too long, overlapping several subwindows in the snap were not recognized in any one of them. In this sense we stamp the false negatives. Some examples of the found meteor trails are shown in Figs. 16, 17 and 18; the numbers are labels of the events. The letter ‘o’ means the original snap segment, ‘d’ means the detected trail (green-highlighted) in this segment. The mark ‘10/1o’ means the original snap of the first part of the meteor trail n. 10 (some meteors leave tracks not only in several subwindows in a snap, but also in two snaps). Examples of long meteor trails that were not recognized at each crisscrossing subwindow are in Figs. 16—3d, 17—10/1d, and 18—15/2d.

Fig. 16

Meteors from January 12/13, 2014, station Ondřejov

Fig. 17

Meteors from station Churáňov

Fig. 18

Meteors from station Churáňov

The survey of the detected events is presented in Table 3 in a chronological order. There are three snap pairs with various trail segments of the same meteor: 1:33:25 and 1:33:55 from December 12, 2014, 3:24:10 and 3:24:40 from March 18, 2015, and 3:33:30 and 3:34:00 from March 25, 2015. Actually, we have 20 different meteors, which can be seen in Figs. 16, 17 and 18.

Table 3 Survey of detected meteors

7.2 False Alarms

Experiments of the first investigated night (12/13 January 2014) have released 383 snaps, which corresponds to the false positive rate of 23%. The second night (23/24 December 2014) experiments put 95 snaps out, which represents the false positive rate of 6.1%. The results of error analysis from both nights are in Table 4, the values in per cents are calculated from the total number of false alarms at the specific night. In the first night, 18 snaps contained double alarms with two different causes, e.g. cloud and airplane, therefore the sum of the false alarms in the first column is 401, but the number of snaps connected with them is 383.

Table 4 The absolute numbers and the percentage of false alarm causes

Some examples of false alarms from these observations caused by a cloud are in Fig. 19, by an airplane in Fig. 20, and by a contrail in Fig. 21. The ratios of the clouds and the contrails are only a rough estimate; in some cases it is practically impossible to distinguish between clouds and contrails, in addition this distinction is for our purpose irrelevant.

Fig. 19

Example of a cloud: a in the original snaps, b found trace

Fig. 20

Example of airplane: a in the original snaps, b found trace

Fig. 21

Example of a contrail: a in the original snaps, b found trace

In one case, moonlight caused a false alarm (see Fig. 22).

Fig. 22

The false alarm caused by moonlight: a in the original snaps, b found trace

Sometimes, the ground objects can cause problems even if they are masked out, see Fig. 23. The chimney looks like a space between dashes from the rotating shutter, or it can prevent subtracting a star, which is hidden in the adjacent snap.

Fig. 23

Ground objects can cause problems even if they are masked out: a, b false space from rotating shutter, c, d star cannot be subtracted, it is behind the ground object in the adjacent snap

False alarms caused by the ground in Table 4 include also interference with the horizon. The combination of clouds, the Earth’s horizon and the camera’s horizon (the circle where the rays from infinity are mapped) and possibly the Moon sometimes creates linear structures that are identified as meteors, see Fig. 24.

Fig. 24

The false alarm caused by moonset behind clouds: a in the original snaps, b found trace

Figure 25 shows another unique case. The hot pixels are so close to one another that they were not erased. They are sufficiently collinear and equally distributed that they have passed all the tests and finally have caused a false alarm.

Fig. 25

Collinear group of hot pixels: a in the original snap, b found trace

An example of a satellite can be seen in Fig. 26. We have two snaps with the satellite image; once (31Dec2015, 20:55:15, camera 0), it was wrongly identified as a meteor, in the other case (31Dec2015, 20:55:45, camera 1), the test correctly rejected it.

Fig. 26

Example of a satellite: a in the original snap, b found trace

The much higher false alarm rate in the first night is caused by two factors. Firstly, the night was cloudy; the clouds were blowing in the wind and were illuminated by the Moon which is the worst case for their elimination. Secondly, the Ondřejov observatory is near Prague and so the airplane traffic and light pollution in this area is much heavier compared to stations in other locations.

8 Analysis of Computing Complexity

The method was originally implemented in Matlab R2015b. We were not satisfied with its computation time, therefore we transcribed it into C++ (Microsoft Visual Studio 2010 Professional). The computation times of individual parts of the algorithm are compared in Table 5.

Table 5 Analysis of computing complexity

The computation was performed on the Intel Core i7-2600K CPU 3.4 GHz processor with 16 GB of memory. The most time-consuming operation is the HT (both with fish-eye correction and without it). The second operation is erasing of hot pixels, here the eccentricity test is much slower than the solitariness test. There are four tests of the lines, one of them is the linearity test. We were not successful in its efficient implementation in the Matlab matrix notation, therefore this test is much slower in Matlab than in C++.

9 Conclusion and Future Work

New equipment and upgrade of the automatic stations require a completely different and novel approach to the processing of observed image sequences. Currently, all data are in digital form and so the approach must also conform to new ways of processing and evaluation. Since vast amount of data obtained from several stations during one night are being processed, it was necessary that we propose a methodology for individual procedures. We tested our partial proposals on real image data and we tried to model various situations that may arise during the search for light traces of meteors. Also, our aim was to carefully sort out the trails originating from a meteor and reduce false alarms to minimum.

There are two remaining elements that could be improved in the algorithms. Firstly, the elimination of clouds is insufficient. In spite of the linearity and high-frequency tests solving this problem, the results are not yet perfect. One possibility for future development is designing additional tests. A similar issue is elimination of airplane trails. Under specific circumstances, it is practically impossible to recognize this trail from the meteor’s trail even for an experienced human observer. Fortunately, these cases are rare and we encountered this situation only once in all the tested sequences. As meteors tend to have more regular dashes, so checking for regularity of the dash positions in trails could be a useful additional test.

The second issue is acceleration of the particular steps. If we need about 11 min per snap, which is captured in half a minute, we ought to distribute the computation among 22 computers (or kernels) to process data on line. The slowest operation is the HT. If we are sure there is just one trail in a subwindow, we can find it by linear regression (or linear with fish-eye correction). Two or more trails in one subwindow would cause their search to fail, but that is very unlikely.

The light curve in a subwindow has a local maximum in time in the snap with a meteor (or other flying object). Much time could be saved if the HT was computed only in the subwindows with such local maximum. The light curve is computed from the mean value of the brightness, which is a very fast operation. The mean value is actually the first-order statistical moment. We also tested variance, skewness, kurtosis and higher moments, but no test had a 100% reliability, therefore we did not use them in a standard fashion. It is possible that another descriptor could be more reliable. Another way to accelerate the process could be a faster hot pixel symmetry test. Still, solutions to related problems depend on the amount of data to be processed. Fortunately, nowadays we are able to handle large data sets using very efficient computing techniques.

We expect applications in other fields that utilize analysis of data sequences, e.g. detection of gamma light flashes, search for particles in cloud and bubble chambers. Other possible application areas include biomedical engineering and organic chemistry materials.