3 Image processing

The purpose of the image processing is to search in the observed field for all the objects which may be detected, in order to identify them and to measure their positions and their fluxes.

Two axes Ox and Oy, where O is the origin of the field, x is along the movement of the strip, and y is perpendicular to x, determine the frame of reference of the observed field. Origin O corresponds to the beginning of the strip in x and to the centre of the strip in y. The columns in y which make up the field are processed one by one in the order of their arrival.

For each column, image processing is achieved in 3 steps:
- extraction of the sky background;
- detection of the objects present in the field and their possible identification by comparison with a reference catalogue;
- measurement of the position and the flux of each detected object.

3.1 Extraction of the sky background

The extraction of the sky background is made either by fitting a polynomial of a degree between 0 and 3 on each column in y, or by using a median filter when the field has strong gradients of flux (e.g. near the Moon, observations of planetary satellites, presence of condensation or hoarfrost, ...). In most cases, there is no significant gradient of flux and the results obtained with both methods are rather similar.

3.1.1 Polynomial fitting

A polynomial is fitted on the flux received on each pixel of the column by least-squares and is substracted from the sky background. Each pixel gives the equation:

$${\rm flux}(i) = a\ [+ bi\ [+ ci^{2}\ [+ di^{3}]]]$$

where i is the number of the pixel in the column, a,b,c,d are the unknowns of the system and brackets indicate optional terms which depend on the degree of the polynomial. Pixels with residuals greater than 3 $\sigma$ (where $\sigma$ is the standard deviation of the residuals) may belong to an object and are therefore removed from the polynomial fit. As changes in the sky background from one column to another are very small, the solution found for one column is used as initial conditions for the next column.

This method is fast but is inconvenient in the case of strong gradients, as mentioned above. However, it can be used in the case of moderate gradients, like those caused by a globular cluster for example.

3.1.2 Median filter

 The principle underlying the classical median filter consists, for each pixel of the image, in searching in a square of $n\times n$ pixels (with $n\geq 3$ and odd) centered on this pixel for the median value of the flux of the pixels, and then substracting this value from the flux of the central pixel. In our case, n must be large enough so that the filter can encompass the largest detected objects. For this reason, a value of n=15 was assumed (Fig. 2).

  

Figure 2: Principle of the median filter

The search for the flux median value of the 225 pixels which make up the filter, when it is done for each of the pixels of the field, requires a lot of computation time and slows down the image processing considerably. That is why a slightly different method was used here. The principle of this method is first to search for the median value Mh(i) of the flux for the pixels along each row i of the filter (i.e. horizontal median), then to search among the 15 values of the flux Mh(i) obtained for the median value M (i.e. vertical median). M is then substracted from the flux of the central pixel of the filter. This method is much faster than the classical method and results from both methods are very close. No pixel rejection is made when searching for the median value of the flux.

3.2 Detection and identification of the objects

After extraction of the sky background, the standard deviation $\sigma$ of the residuals of the pixels of the processed column is calculated. Pixels with residuals greater than 3 $\sigma$ may belong to an object. In order to avoid as far as possible identifying some possible parasitic effects with actual objects, an object is defined by at least 2 consecutive pixels with residuals greater than 3 $\sigma$. Detection threshold may be lowered down to 2.5 $\sigma$ when searching for faint objects. Limits of the object are searched for in the processed column and in the following ones, then a detection window is defined surrounding the object. A security margin of 2 pixels in each direction was added between the edges of the object and the limits of the window, in order to ensure that the whole object is contained in the window. True local sidereal time (equal to the apparent right ascension when the object is on the meridian) and an approximation of the apparent declination of the centre of the window are known and compared with the apparent positions of the stars in a reference catalogue. If there is agreement to within $\pm$3 pixels, the object is identified with the star giving the best agreement. If not, a temporary identification is given.

3.3 Position and profile of the objects

The search for the photocentre of the object is made by using a two-dimensional Gaussian flux distribution:

$$\begin{eqnarray} \phi(x,y) & = & \frac{\Phi}{2\pi\sigma_{x}\sigma_{y}\sqrt{1-\rh... ...a_{x}}\right) \left(\frac{y-\sigma_{y}}{\sigma_{y}}\right)\right]}\end{eqnarray}$$ (1)

where $\Phi$ is the total flux of the object, $(\mu_{x},\mu_{y})$ is the centre of the Gaussian, $\sigma_{x}$ and $\sigma_{y}$ are the standard deviations in x and y of the Gaussian and $\rho$ is the correlation coefficient between the two dimensions of the Gaussian. Equation (1), when written for each pixel of the object, leads to a non linear system which has to be solved. That is why for each pixel (i,j) a linearization is made:

$$\begin{eqnarray} \phi_0(i,j)-\phi(i,j) & = & \frac{\partial\phi}{\partial\Phi}\D... ..._{y}}\Delta\sigma_{y} +\frac{\partial\phi}{\partial\rho}\Delta\rho\end{eqnarray}$$ (2)

where $\phi_{0}(i,j)$ is the flux received by pixel (i,j), $\phi(i,j)$ is the corresponding flux calculated by formula (1), and $\Delta\Phi,\ \Delta\mu_{x},\ \Delta\mu_{y}, \ \Delta\sigma_{x},\ \Delta\sigma_{y},\ \Delta\rho$ are the unknowns of the system. This system is solved by least-squares and by successive approximations. This means that the barycentre of the image is first calculated and then used as initial conditions for fitting a circular Gaussian (i.e. $\sigma_{x}=\sigma_{y}$ and $\rho=0$). If the circular Gaussian converges, fit of a non-circular Gaussian is attempted. If the non-circular Gaussian does not converge, the assumed flux distribution is the one given by the circular Gaussian. If the circular Gaussian did not converge, the photocentre is assumed to be the barycentre of the image.

  

Figure 3: Undersampling case: surfaces S1 and S2 are clearly different. So here, the difference between the flux received on the pixel (surfaces S+S1) and the corresponding flux given by the Gaussian (surfaces S+S2) is certainly significant

In most cases, the system is undersampled, that is to say, the size of the pixel (1.65 or 1.51'') is too large as compared with the standard deviation of the Gaussian (about 0.7 pixel in good conditions). In these cases, the difference between the received flux on a pixel and the corresponding flux given by the Gaussian for this pixel is significant (Fig. 3). This phenomenon leads to a bias in the determination of the photocentre of the star which depends on the position of the photocentre with respect to the pixel. To solve this problem, a more exacting calculation was made, by replacing in formula (2) $\phi(i,j)$ with $\phi_{\rm cor}(i,j)$, with:

$$\phi_{\rm cor}(i,j) = \int_{i-\frac{1}{2}}^{i+\frac{1}{2}} \int_{j-\frac{1}{2}}^{j+\frac{1}{2}}\phi(x,y)\;{\rm d}x{\rm d}y$$

wherein $\phi(x,y)$ is given by formula (1). In order to simplify the calculation of $\phi_{\rm cor}$, a Taylor series expansion of $\phi(x,y)$ to the 2$^{\rm nd}$ order in x and y was made in the vicinity of (i,j), with the condition {$-1/2\le x-i\le 1/2$, $-1/2\le y-j \le 1/2$}. In this expansion, terms depending on $\rho$ are considered to be negligible. The result is:

$$\phi_{\rm cor}(i,j)=\phi(i,j)\left\{1+\frac{r_{x}^{2}-1}{24\sigma_{x}^{2}} +\frac{r_{y}^{2}-1}{24\sigma_{y}^{2}}\right\}$$

with:

$$r_{x}=\frac{i-\mu_{x}}{\sigma_{x}},\ r_{y}=\frac{i-\mu_{y}}{\sigma_{y}}$$

which is equivalent to applying a correction to $\phi(i,j)$.

Another solution to the problem of undersampling, which may be combined with the preceding solution, is to defocus star images, leading to a larger FWHM (full width at half maximum) of the Gaussian, as Stone et al. (1996) did on the Flagstaff Astrometric Scanning Transit Telescope.

Dynamical range for the CCD receptor is about 7 mag. For bright objects ($V\ \le\$8.5 in most cases), some pixels of the object are saturated. These pixels are not used for fitting the Gaussian, in order to avoid as far as possible any loss in precision on the determination of the photocentre and any bias on the measurement of the flux of these objects.

On the other hand, if multiple objects have been detected (i.e. several maxima appear in the square of pixels surrounding the main image), several Gaussians are fitted simultaneously in order to enable the position and the flux of each object to be measured independently, when separation between the objects is larger than 5''.

3.4 Case of high declinations

In drift scan mode, the rate of motion of electric charges in the detector is matched with the mean rate of transit of the star image in the field. The differential transit rate of the star image between the top and the bottom of the CCD detector increases with declination. In the case of large declinations (positive or negative), some limitations appear because the differential transit rate between the centre and the borders of the field is too high. Another problem is the fact that the star in front of the detector follows a curved path (due to the curvature of the parallel, which increases with declination). In the case of large declinations, both problems result in distorsion of the star's image and loss of precision in its measurement, which is greater if the star is near the borders of the field. Accordingly, stars for which the elongation of their image along the x axis is greater than $6\;\sigma$ (where $\sigma$ is the standard deviation of the star measurement) cannot be used for reduction and are rejected. With the 1024$\times$1024 detector, the problem becomes significant for declinations above about $45^{\circ}$. At $\delta=60^{\circ}$, about 50% of the field is rejected. At Valinhos, this problem is less critical because the size of the detector is twice as small. In all cases, declinations above $\pm 77.8^{\circ}$ are out of reach because of electronic limitations.