Measurement uncertainty of the computed center of gravity from pixel image

Question

Establishing the position of celestial bodies (stars, and other objects) is based on imaging technology. Depending on the resolution of the imaging sensor used, an object is identified by a group of neighboring (bright) pixels. These pixels surround the center of gravity of the detected celestial body.

In image processing, the central position of a round object (a sphere) is generally determined by computing the center of gravity of the pixels surrounding the detected object. The weights are the intensities of the pixels in the vicinity of the brightest (central) pixel.

My question: how is the variance of this center-of-gravity estimate being calculated, within astronomy?

Answer 1

Every pixel value $S_i$ on the detector at $\vec x_i$ has some error $N_i$: CCDs for example have a background noise $N_\text{bkg}$ from read-out electronics, thermal noise, and sky background, plus a Poissonian photon noise $N_S=\sqrt{S}$. In many cases this noise follows a Gaussian distribution reasonably well. After subtracting the background, a position measurement

$$ \vec x = \frac{\sum{S_i \vec x_i}}{\sum{S_i}} $$

has the uncertainty

$$ \operatorname{std} \vec x = \left[\sum_i \left(\frac{\partial x}{\partial S_i} N_i\right)^2 \right]^{1/2} = \frac{\left[\sum_i \left( \sum_j\frac{(x_i - x_j) S_j}{\sum_k S_k} N_i \right)^2 \right]^{1/2}}{\sum_k S_k} $$

$$ N_i^2 = N_\text{bkg}^2 + N_{S_i}^2 $$

using Gaussian error propagation ... I hope I have run the math correctly. The way I've written it is a bit weird, but shows one important property of this method: If you look at the sum over $j$, you can see how the noise of a pixel is basically weighted by the pixel's distance to the centroid result. The same error on a pixel value has more impact if the pixel is farther away from the centroid.

Better methods account for the uncertainty $N_i$ on pixel values in the first equation already. You can do this by introducing additional weights in your centroid or go to fitting models, which is the “usual” method in astrometry, as far as I have seen.

This more involved approach to measuring positions uses the point spread function $P(\vec x_i - \vec x_\text{obj})$ of the instrument under given observing conditions. The model can be an approximation, e.g. a Moffat function, or an empirical model build from unconfused bright stars in the image, e.g. a spline interpolation. For point sources, the typical least-squares fit of position and flux of the model to an image easily produces results close to the statistical optimum of parameter estimation in terms of uncertainty. With our modern computing power the easiest way to get the uncertainty for a given model and data noise is often a bootstrapping algorithm.

Of course extended objects require a bit more work in the model, for example assumptions on their shape, as you have given in your question.

Answer 2

You might take a look at this 2006 paper by Thomas et al. on centroiding algorithms for astronomical adaptive-optics (AO) systems, which includes a detailed discussion of error estimates for the centroid position using different algorithms. The approach you describe in your answer corresponds to what they call "simple centroid" (Section 3); they refer to a book chapter by Rousset (1999) for the detailed analysis (which I believe includes Poisson and readout-noise contributions, and so is not identical to your result).

More generally, the "center-of-gravity" approach seems to be used in situations where one needs fast, computationally cheap estimates, such as in AO systems (where a stellar centroid needs to be determined many times a second), or as a crude first guess to provide a starting point for a more sophisticated analysis. Post-observing analysis of astronomical images generally uses more complex/sophisticated approaches, depending on whether the source is a star or other point sources versus an extended source, whether you have an accurate model of the point spread function (which may be non-circular), deblending of neighboring sources, what the noise characteristics of your data are, etc.

In practice, I'd guess that most such analyses use some kind of nonlinear least-squares or maximum-likelihood analysis that involves fitting a model to the data. The errors in the fitted model parameters (including the centroid position) can be derived from simplistic assumptions about the fit landscape (e.g., Levenberg-Marquardt and other gradient-based minimization algorithms sometimes provide covariance matrices based on treating the local $\chi^{2}$ landscape as a parabola), from bootstrap resampling, or from Markov Chain Monte Carlo approaches. This may be supplemented by running simulations of the fitting process on artificial images of simple models of stars or galaxies, to get some quasi-empirical estimates or corrections for the centroid (and other parameter) uncertainties.

Answer 3

Dunno about "...generally determined..." but in similar work I've performed (and witnessed) a 2-dimensional spline fit is performed on the pixel intensity data to determine the power centroid to sub-pixel resolution. As you noted, we make the reasonable assumption that the object is close to spherical and close to having azimuthal constant density (that is, density may vary with radius but not angle).

The uncertainty (variance) in this calculation is usually calculated by applying standard statistical methods to the observed variances in the received signal in each pixel, after accounting for jitter in the line-of-sight position. (and of course accounting for electronics noise, etc). Essentially, one can look at the variance in the spline-calculated peak over N frames, or look at the variance in all the pixels, determine their contribution to the spline fit, and weight their contributions accordingly.

Answer 4

This interesting problem is clearly relevant in astronomy as well. Within medical imaging, determining the center of an object also comes into play in several, special application areas.

I posed the question here to become informed as to which formula's are being used within astronomy: How is the uncertainty of the 'COG' being estimated, in the presence of (additive) measurement noise?

In 2002, we derived a general formula for the variance of the estimated center of gravity, in the presence of normally distributed additive noise, associated with each pixel value. The reference is as follows:

H.C. van Assen, M. Egmont-Petersen, J.H.C. Reiber. "Accurate object localization in gray level images using the center of gravity measure: accuracy versus precision", IEEE Transactions on Image Processing, Vol. 11, No. 12, pp. 1379-1384, 2002.

I will first give the general formula, which only assumes that the center of the window surrounding the object is $(0,0)$. The object does not have to be positioned centrally for this formula to hold.

Define the additive measurement noise associated with each pixel as $\epsilon \sim U(0,\sigma^2)$, with $\sigma$ being its standard deviation.

Define a squared (sub)image ${\cal W}$ of dimensions $(d+1) \times (d+1)$ pixels, with coordinates, $x=-\frac{d}{2},\frac{d}{2}+1,\ldots,\frac{d}{2}$, and $y=-\frac{d}{2},\frac{d}{2}+1,\ldots,\frac{d}{2}$. Let ${\bf w}_{x,y}$ be the true pixel intensity of pixel $(x,y)$ in absense of any noise, in ${\cal W}$. Define the signal-plus-noise image $W$ as: $w_{x,y}={\bf w}_{x,y} + \epsilon$. Define the total number of pixels in the (sub)image as $N=(d+1)^2$, which includes the central 0-th row and 0-th column of $W$.

The values of $w_{x,y}$, with ($w_{x,y} \geq 0$), are the ones actually being observed. Close to the center $(x=0,y=0)$ a bright object of interest has been located.

The estimated center of gravity $cog$ of this object is computed as follows: $$ \widehat{cog}(x,y)=\left(\dfrac{\sum_{x,y} \; x\,w_{x,y}}{\sum_{x,y} \; w_{x,y}},\dfrac{\sum_{x,y} \; y\,w_{x,y}}{\sum_{x,y} \; w_{x,y}}\right) $$ where $x$ and $y$ are running indices such that each pixel in $W$ enters the calculation of a sum (sigma) exactly one time. $cog$ is the measurement that is influenced by measurement noise.

Using the delta rule two successive times, we derived a general approximate formula for the variance of the cog, given a known noise level.

Define $x2$ as: $$ x2 = \sum_{x} \; \sum_{y} \; x^2 $$ and similarly $y2$ as: $$ y2 = \sum_{x} \; \sum_{y} \; y^2 $$ Finally, the average 'weight' (average pixel intensity) is given by: $$ \hat{\mu}_w = N^{-1} \; \sum_{x} \; \sum_{y} \; w_{x,y} $$ Let $\widehat{cog}(x)$ denote the estimated center of gravity $x$-coordinate and $\widehat{cog}(y)$ the estimated center of gravity $y$-coordinate.

The MLE derived variance estimates of $\widehat{cog}(x)$ and $\widehat{cog}(y)$ are: $$ \text{var}(\widehat{cog}(x))=\left(\frac{\sigma^2 \; x2}{\left[N \; \hat{\mu}_w \; \widehat{cog}(x)\right]^2} + \frac{\sigma^2}{N \; (\hat{\mu}_w)^2} \right) \; (\widehat{cog}(x))^2 $$ and $$ \text{var}(\widehat{cog}(y))=\left(\frac{\sigma^2 \; y2}{\left[N \; \hat{\mu}_w \; \widehat{cog}(y)\right]^2} + \frac{\sigma^2}{N \; (\hat{\mu}_w)^2} \right) \; (\widehat{cog}(y))^2 $$ It turns out that when the true cog is precisely $(0,0)$, then the limit (simplified) formulas for the cog-variance hold: $$ \lim_{cog(x) \to 0} \; \; \lim_{cog(y) \to 0} \; \text{var}(\widehat{cog}(x)) = \frac{\sigma^2 \; x2}{(N \; \hat{\mu}_w)^2} $$ and $$ \lim_{cog(x) \to 0} \; \; \lim_{cog(y) \to 0} \; \text{var}(\widehat{cog}(y)) = \frac{\sigma^2 \; y2}{(N \; \hat{\mu}_w)^2} $$ as the second (additive) variance term within the outer parentheses then vanishes.

My recent simulations show that above $95\%$ of the actual variance results from out defined variance formula, as presented here. This simulation result also holds when the true cog deviates more than one coordinate from the central position $(0,0)$.

Simulation graphic showing the accuracy of the variance estimation will be added to this answer, one of the coming days.

Multiplicative noise

Poisson noise occurs in CCD-cameras, its influence can especially be observed in gradient transitions of image intensity from very dark areas to really bright ones. The magnitude of Poisson noise is known to be proportional to signal intensity. In this case, multiplicative noise is present in the observed image: $$ w_{x,y}={\bf w}_{x,y} \, \, \cdot \alpha \, \cdot \, \epsilon $$ with $\alpha$ being a proportionality constant and $\epsilon$ the normally distributed noise term (which approximates a Poisson distribution well for moderate $w_{i,j} >0$).

Performing a simple $\ln(\cdot)$ transformation, $lw_{i,j}=\ln(w_{i,j})$ yields a transformed image $lw$ perturbed by additive noise. Subsequently, the cog-variance estimator can be applied to the $lw$-image.