3 Technical aspects

3.1 Bad pixels

StarFinder includes an automatic procedure to repair known bad pixels, which are replaced with the median of the good data points in a suitable neighborhood. The use of this procedure is particularly recommended in the extraction of the PSF. Indeed the stellar images, which will be combined to form the PSF, are centered with sub-pixel accuracy by means of an interpolation technique and the presence of bad pixels may produce unpredictable interpolation errors. Even if replaced, the bad pixels are excluded from any further computation concerning the correlation coefficient and the fitting process.

3.2 Background estimation

A reliable estimate of the background is necessary to define the detection condition and to compute the correlation coefficient of the presumed stars with the PSF.

A straightforward technique to estimate the background is to smooth the image by median filtering, replacing each pixel with the median computed over a suitable neighborhood, of size comparable to the characteristic width of a stellar image. This method tends to over-estimate the background underneath strong peaks. A more accurate approximation, even in crowded fields, is described in Bertin & Arnouts (1996). The image is partitioned in sub-regions arranged in a regular grid and a local estimate is calculated for each sub-patch by means of an IDL implementation of the DAOPHOT SKY routine (Stetson 1987), due to Landsman (1995). This array of sky measurements is smoothed by median filtering and interpolated onto the same grid of the input image.

It should be stressed that the background computation is unavoidably affected by the presence of bright sources. In general a more accurate estimate can be obtained after the analysis of the stellar field, when the most contaminating sources are known and can be subtracted.

3.3 Noise estimation

The estimate of the noise is useful to define the detection threshold and to compute the formal errors on the retrieved astrometry and photometry.

The overall effect of photon and instrumental noise can be computed if the required parameters (detector gain, read-out-noise, dark current, etc.) are known. Otherwise an estimate of the mean background noise can still be obtained by means of histogram fitting techniques (Almoznino et al. 1993; Bijaoui 1980). Assuming that the intensity of the sky radiation is distributed normally around a typical value, the histogram of the observed intensity levels should be quite similar to a Gaussian distribution, whose mode and standard deviation represent respectively the sky value and the associated noise. Actually the contamination due to the stellar sources produces a high-intensity tail and an artificial broadening of the histogram, which prevents an accurate estimate of the background noise. This problem can be partially overcome by removing the signal from the image, leaving only the pixel-to-pixel variations associated to pure noise: a reasonable estimate of the signal to subtract for this purpose can be obtained by smoothing the data with a median filtering technique. After this operation the histogram is symmetric around its mode and the background noise standard deviation can be estimated by means of a Gaussian fit to the histogram itself.

3.4 Correlation coefficient

False detections, associated to noise spikes or residual PSF features of bright stars, are recognized on the basis of their low correlation coefficient with the PSF, which represents a template for each true star in the image. The correlation coefficient (Gonzalez & Woods 1992) is defined as

$$% c(a,b) = \frac{\sum_{x,y}\left[ i\left( x,y\right) -\bar{\imath... ...ht]^{2}}\sqrt{\sum_{x,y} \left[ p\left( x-a,y- b\right) - \bar{p}\right] ^{2}}}$$ (2)

where i(x,y) and p(x,y) are the object and the PSF respectively, $\bar{\imath}$ and $\bar{p}$ are the corresponding mean values. Maximizing the correlation coefficient as a function of the offset (a,b) yields an objective measure of similarity. After maximizing c(a,b) for integral offsets, it is possible to repeat the procedure for sub-pixel shifts, improving the positioning accuracy.

The correlation coefficient is computed on the core of the PSF: typically the central spike of the diffraction pattern is considered, out to the first dark ring. A fair correlation threshold must be fixed in order to discriminate and reject unlikely detections, without losing faint stars contaminated by the background noise; a value of 0.7 or 0.8 is acceptable in most cases.

The correlation coefficient represents also an effective tool to select the stars with the highest photometric reliability, since generally a very high correlation value is associated to resolved single sources.

3.5 Saturated stars

Saturated stars provide precious information on the PSF halo, so it may be useful to include them, appropriately reconstructed, in the PSF extraction process. In addition, the repaired saturated stars can be recognized by the star detection algorithm and their contribution be taken into account during the analysis of near and fainter sources.

The core of a saturated star is replaced with a shifted scaled replica of a preliminary estimate of the PSF. The repaired star is defined as

$$% i_{\rm rep.}(x,y)= \left\{ \begin{array}{ll} i(x,y) & \mbo... ...m] f p(x-x_{0},y-y_{0}) & \mbox{otherwise} \end{array}\right.$$ (3)

where T is the upper linearity threshold of the detector. The position of the center (x₀,y₀) is estimated by means of the maximization of the correlation coefficient, which is not sensitive to the intensity levels. The scaling factor f is calculated with a least squares fit to the wings of the saturated star. Of course the saturated pixels are excluded from all the computations. The background is temporarily subtracted in order to prevent affecting the repair process.

This procedure has been applied to the brightest star (IRS 7) of the Galactic Center image shown in Sect. 4.2. This source is not saturated in the original data, but it has been artificially corrupted by an upper cut at half maximum. The repair procedure is able to reconstruct the original peak with an error of $\sim$5%, imposing a positioning accuracy of 1/2 pixel.

3.6 Fitting procedure

A sub-image centered on the star of interest is extracted and approximated with the model

$$% h(x,y) = s_{0}(x,y)+\sum_{n=1}^{N_{\rm s}}f_{n}p\left( x-x_{n},y-y_{n}\right) +b_{0}+b_{1}x+b_{2}y$$ (4)

where s₀(x,y) is the fixed contribution of known stars outside the sub-image support, $N_{\rm s}$ is the number of point sources within the sub-image, x_n, y_n, f_n are the position and flux of the n-th source, p(x,y) is the PSF and b₀, b₁, b₂ are the coefficients of a slanting plane representing the local background. It should be stressed that the retrieved astrometry and photometry are referred to the absolute centering and normalization of the PSF array. The optimization of the parameters is performed by minimizing the least squares error between the data and the model. If the noise is known, it is possible to weigh the data by their inverse standard deviation and obtain a statistically optimal fit (Beck & Arnold 1977). In this case, formal error estimates on the parameters can be obtained (Bevington & Robinson 1994). The optimization is performed by means of an iterative Newton-Gauss technique with linearized Hessian (Beck & Arnold 1977; Luenberger 1984). The Moore-Penrose generalized inverse concept (Rust & Burrus 1972; Bendinelli et al. 1987; Lorenzutta 1967) is applied to invert the Hessian matrix at every step. The iterations are stopped when the parameters approach a stable value. The major difficulty of a Newton-like method is represented by the computation of the model derivatives with respect to the parameters, some of which (stellar positions) yield non-linear dependence. For this purpose the Fourier shift theorem is applied:

$$% p(x-x_{n},y-y_{n}) = {\rm FT}^{-1} \left[{\rm FT} \left [p(x,y) \right]{\rm e}^{-i2\pi (ux_{n}+vy_{n})/N}\right]$$ (5)

where FT represents the discrete Fourier transform operation, N is the sub-image size and u, v are spatial frequencies. The derivative with respect to x_n is then

$$% { \frac{\partial p\left( x-x_{n},y-y_{n}\right) }{\partial x_{n}} = }$$

$${\rm FT}^{-1}\left[ {\rm FT}\left[ p\left( x,y\right) \right]{\rm... ...-i2\pi \left(ux_{n}+vy_{n}\right) /N}\left(-i \frac{2\pi u}{N} \right)\right] =$$

$${\rm FT}^{-1}\left[ \left( -i\frac{2\pi u}{N}\right) {\rm FT}\left[ p\left(x-x_{n},y- y_{n}\right) \right] \right]$$ (6)

and requires in practice an interpolation of the PSF to compute p(x-x_n,y-y_n). In principle this interpolation-based method can only be applied to Nyquist-sampled images and this is currently the main limit of the algorithm. A similar technique has been described by Véran & Rigaut (1998).

3.7 Sampling and interpolation

Figure 4: Synthetic field image with 1000 sources. The PSF Strehl ratio is $\sim$40%. The display stretch is logarithmic

A band-limited function in one dimension is Nyquist-sampled if the step size fulfills the condition

$$% \Delta x\leq \frac{1}{2}f_{\rm c}$$ (7)

where $f_{\rm c}$ is the so-called cut-off frequency of the spectrum. In this case it is possible to reconstruct the original continuous function from a set of equally-spaced samples by means of the so-called sinc interpolation (Mariotti 1988). For a nearly diffraction-limited PSF, the critical sampling condition is generally stated by saying that its FWHM must contain at least two pixels. Many accurate and efficient interpolation schemes exist to perform the PSF shift required by the fitting procedure. Probably the most simple is a straightforward application of the Fourier shift theorem, but one may also use bicubic splines or the fast sinc interpolation algorithm described in Yaroslavksy (1997). In StarFinder we have applied the cubic convolution interpolation method (Park & Schowengerdt 1983), which approaches very closely the optimal sinc interpolation for Nyquist-sampled data; this algorithm is implemented in the IDL function INTERPOLATE.

Even though the interpolation of the PSF array is allowed only on Nyquist-sampled data, the cubic convolution technique seems to produce acceptable results even with marginally under-sampled images. Tests performed on Airy diffraction patterns, with a sampling step twice as large as the critical sampling step size, indicate that the interpolation-induced oscillations amount to a few percent of the image peak, as opposed to 10-20% of other interpolation techniques like Fourier shift or bicubic splines.