Star detection

Star detection is an algorithm that finds stellar objects on a CCD frame. The input of the algorithm is a calibrated CCD frame F; that is, the bias, dark and flat corrections have already been applied to the source image data. The output of the algorithm is a table of x and y coordinates that correspond to the centers of detected stellar objects. For practical use, the robustness of the detection is a crucial property, because the signal to noise ratio of real images tend to be low, and there are often artifacts present.

The algorithms that are in the C-Munipack software follow the algorithms from the Stetson’s FIND routine - see [stetson11] for the documentation of the original code.

Estimation of overall background level and noise level

In the prologue of the star detection, the overall background mean level and noise level is determined. These values are required in the following phases. The robust mean algorithm, described in the section Robust mean, is used to determine the mean level and standard deviation of the background. It is assumed, that majority of pixels are not stars; in this case the pixels that are covered by stars are rejected by the robust mean as outliers, but if the field was almost covered by stars, the background estimate would be biased.

Not all pixels are included in the computation of the robust mean, the program filters out pixels that belongs to the user-defined border and also bad and overexposed pixels.

Fitting Gaussian function

Because of the noise, the shapes of stars are not regular. Especially for faint stars, the maximum pixel need not be located at its center. When the telescope was slightly de-focused during the exposure, the profiles of stars has a depression in its center caused by a projection of the secondary mirror on the detector. Insufficient filtering at this stage causes multiple detections of a single object.

The low-pass filter is implemented by fitting a function y to each pixel using several pixels in its neighborhood. The function y is a linear combination of a sampled two-dimensional Gaussian function g(i,j) that represents an image of a stellar object and a constant function that represent a local background level. The linear least squares method is used to determine the two coefficients h and s of fitted function y:

(1)y(i,j) = h\,g(i,j)+s = h\,e^{-\frac{(i-i_0)^2+(j-j_0)^2}{2c^2}}+s

where i_0 and j_0 are coordinates of a central pixel and c is computed from the configurable parameter FWHM, which determines the width of a Gaussian function and the filter’s frequency response. The parameter c is related to the FWHM parameter according to:

(2)FWHM = {2\sqrt{2 \log 2}}\,c

The linear least squares find values of free parameters h and s which minimizes a sum of squared residuals between original image data F(i,j) and its model y(i,j)

(3)S(h,s) = \sum{(y(i,j)-F(i,j))^2} \rightarrow min.

The Gaussian function as well as the constant function are contiguous and indefinite, for practical implementation we have to restrict the integration space to the pixel’s neighborhood. We determine the half-length of a filter using the FWHM parameter:

(4)\text{filter half-length} = \max(2, \lfloor 0.637 * \text{FWHM} \rfloor)

For practical use, we define n as a width and height of a square matrix as n = 2 * \text{filter half-length} + 1 and also a value of N as a number of pixels in the pixel’s neighborhood N = n^2.

From a solution of a linear least squares we get the following formula for the coefficient h that multiples the Gaussian function.

(5)h = \frac{N\,\sum{gf}-\sum{f}\sum{g}}{N\,\sum{g^2}-(\sum{g})^2}

Equation (5) is applied to each pixel of a source frame F independently, leaving out pixels that are too close to the border, less that filter half-length. Computed values of the h coefficient for each pixel are stored in a new auxiliary frame H.

Finding star candidates

The object candidates are found by searching for local maxima in the auxiliary frame H. The neighborhood of each pixel at coordinates x and y is checked out to the radius equal to the filter half-length. If the pixel in the center has a greater value than all pixels in its neighborhood, we have found a candidate.

The candidates are subjected to several subsequent tests. Each test was designed to be sensitive to specific class of artifacts.

Checking minimum brightness

The height of the fitted two-dimensional Gaussian function stored in the filtered frame H of a candidate with local maximum at coordinates x and y must be equal to or greater than minimum height, hmin:

(6)hmin = \text{THRESHOLD} * relerr * \sigma_{F}^2

where THRESHOLD is the detection threshold parameter, adjustable in the configuration of the program.

The relerr term in the equation is determined as:

(7)\frac{1}{relerr} = \sqrt{\sum{G^2} - \frac{(\sum{G})^2}{N}}

where N is a number of pixels in the matrix G.

The \sigma_{F}^2 term is an estimation of random noise per pixel in ADU. The random noise has two independent sources. The first one is the Poisson statistics that applies to a number of photons that hit the detector during the exposure (number of discreet events that occurred in a fixed time). The variance \sigma^2 of the Poisson random variable with mean value of \lambda is \sigma^2 = \lambda. If we know the number of photons per ADU, the gain p of the detector, and the mean level s of the background, the mean number of photons \lambda that hit the detector is p.s.

The variance of the random error per pixel in ADU \sigma_{ADU}^2 is determined from the formula

(8)\sigma_{ADU}^2 = \frac{s}{p}

The second source of noise is the preamplifier and A/D converter \sigma_{RN}^2, usually referred to as the read-out noise. This parameter should be given by a manufacturer in the camera’s data sheet, however, this value is stated for a single frame and when a processed frame has been made by summing or averaging of multiple raw frames, a corresponding transformation of the specified value of the read-out noise must take place. This applies namely when one is processing DSLR images or combined CCD images.

If the processed frame has been made as a sum of N raw frames, then

(9)\sigma_{RN}^2 = (\text{READOUT\_NOISE})^2 * N

where READOUT NOISE is a read-out noise value in A/D units per frame (from the data sheet).

If the processed frame has been made as a average (arithmetic mean) of M raw frames, then

(10)\sigma_{RN}^2 = (\text{READOUT\_NOISE})^2 / M

These sources are independent, therefore from the error propagation rules we determine the overall noise \sigma_{F}^2 as:

(11)\sigma_{F}^2 = \sigma_{ADU}^2 + \sigma_{RN}^2

Sharpness of an object

This test is intended to leave out diffuse objects, traces of cosmic particles, etc. The sharpness value of a candidate must be within configurable limits. Sharpness is defined as a ratio of height of best-fitting Gaussian function and best-fitting sampled two-dimensional \delta function.

(12)\delta(i,j) = \begin{cases} 1, & \textrm{if}\,i=i_0\,\textrm{and}\,j=j_0 \\ 0, & \textrm{otherwise} \end{cases}

The first value, the height of best-fitting Gaussian function has already been determined in advance – it is stored in the auxiliary frame H. The procedure of fitting the delta function is the same as for the two-dimensional Gaussian function. We model the image as a linear combination of the delta function \delta and the constant function.

(13)S(d,t) = \sum{(d.\delta(i,j)+t-F(i,j))^2} \rightarrow min.

We limit the fitting of the delta function to the same number of neighboring pixels as we did for fitting the Gaussian function. By solving the linear equations to find a value of free parameters d and t that minimizes the sum of squared residuals between the model y(i,j) and original frame F we get a height of best-fitting delta function as a value of d:

(14)d = \frac{N\,\sum{\delta f}-\sum{f}\sum{\delta}}{N\,\sum{\delta^2}-(\sum{\delta})^2} = \frac{N\,f(i_0,j_0)-\sum{f}}{N-1} = f(i_0,j_0) - \frac{\sum{f}-f(i_0,j_0)}{N-1}

The last form of the equation is easy to implement. It says: take the central pixel and subtract the average of all surrounding pixels.

Roundness of an object

This test rules out candidates that have a width along the x-axis much different than along the y-axis. Such local maxima do not correspond to real stars, but to artifacts caused during the read-out of the CCD chip.

The roundness parameter is defined as the relative difference between height of the best-fitting Gaussian in x and y axes. Please note, that unlike the previous steps, we fit the image data to a one-dimensional Gaussian function. Therefore, we need to reduce the dimensionality of the image data. The distribution of the signal is modeled after a 2D Gaussian function. We take advantage of a fact that an integral of the 2D Gaussian function G(x,y) in x is 1D Gaussian function G(y) and vice versa.

First, we sum the image data in the neighborhood of the local maxima by rows to a vertical profile of an object and we fit the result with a linear combination of a one-dimensional Gaussian function and a constant function. From the least squares fit we get the height of Gaussian h_x. Second, we sum the image data by columns to get a horizontal profile of an object and by the same way we get a value h_y. The roundness parameter is derived as the difference between those two values divided by their average.

(15)\textrm{ROUNDNESS} = 2\,\frac{h_x - h_y}{h_x + h_y}

If the heights of both Gaussian functions are equal, the value of roundness is zero. If one the Gaussian functions has zero height, the value is 2 or -2. The default allowed range for the roundness is -1 and 1.

Finding center of the object

Up to this point, we have worked with the data as if the center of the object was located in the center of the brightest pixel. In this step, we derive the position of the object by estimating an offset vector of the centroid of the brightness enhancement to the brightest pixel.

The offset vector is derived in two subsequent steps, its d_x and d_y components are computed separately. The value of d_x is computed by fitting the first derivative of a one-dimensional Gaussian function G'(x) to the difference between image data F(x) and the best-fitting Gaussian function G(x). We use the same horizontal profile of the object as we used in the computation of the roundness parameter and we subtract the best-fitting Gaussian function h_xG(x) that we have already derived. Then we fit the result of subtraction by the linear combination of the function G'(x) and a constant function using the linear least squares method.

(16)S_x(j_x,k_x) = \sum{((j_x.G'(x)+k_x)-(F(x)-h_x.G(x)))^2} \rightarrow min.

We find the value of j that minimized the sum in equation (16). The value is then used to estimate the offset between the central pixel and the centroid of the brightness enhancement using the equation (17). Proof?

(17)d_x = \frac{j}{1 + |j|}

The same method is used to derive the offset d_y in the vertical direction.

The object’s center \vec{o} is computed as a sum of coordinates of central pixel \vec{x} and the offset vector \vec{d} = (d_x,d_y).