- This topic has 13 replies, 4 voices, and was last updated 4 years, 2 months ago by
.
-
Posts
-
I am trying to resolve some problems related to calculating the SNR of a photometric measurement when all you have is a single calibrated image, i.e. no way of extracting the individual noise components. The reason I am doing this is that Astrometrica gets this wrong. In fact it reports two different SNRs in two different places and neither corresponds to the SNR I get from AstroImageJ.
The basic approach seems to be to measure the total star count in an inner circular aperture and the sky mean and standard deviation in a concentric outer annulus. The sky standard deviation then represents all of the noise components except for the star’s photon noise. There is then a formula for calculating SNR in one of the documents for VPHOT (Page 5 of here) but I am pretty sure that this is wrong and the brackets are in the wrong place. My guess is that it should be something more like:
SNR = S / sqrt( (S/G) + N*sigma^2 * (1+1/M) )
where S is the star count, sigma is the sky standard deviation, N is the number of pixels in the star aperture, M is the number if pixels in the sky annulus and G is the camera gain (ADUs/photon). I have tried this and it kind of works but the numbers are different to AIJ and Astrometrica so I now have four different SNR values! Also, if I stack images, sigma goes down as sqrt(N) where N is the number of images but the SNR does not go up as sqrt(N). Why not?
This estimator must have been implemented loads of times. Does anyone have a definitive reference which covers how to do this?
I’m getting very frustrated with this and would be pulling my hair out if I had any…
Nick,
You’ve got me attempting to do undergrad statistics at 2am after a couple of beers, so this may be a pile of nonsense. But the formula I get is as follows…
* By definition, SNR = Signal / sqrt(Sky_Noise^2 + Photon_Noise^2)
* Define the sky noise in an individual pixel to be sigma_sky. If the standard deviation of M sky pixels is measured to be sigma, then your best estimate of the sigma_sky = sigma*sqrt(M/(M-1)). The sqrt() term comes about because sigma is a slight underestimate of sigma_sky, since the same sample of M sky pixels was also used to calculate the mean sky level relative to which their standard deviation was measured. For large M, sigma_sky=sigma.
* The total noise in the summed brightness of all N pixels within your aperture is sigma_sky*sqrt(N) = sigma*sqrt(NM/(M-1)).
* This equals Sky_Noise in the SNR equation above.
* Now for photon noise. Define the number of photons collected from the star to be N_photon. By the definition of the gain, N_photon = S/G. The noise in N_photon is sqrt(S/G) since it is a Poisson process. Thus, the noise in S is sqrt(S/G)*G = sqrt(SG).
* Putting it all together…
SNR = S / sqrt(SG + sigma^2 * NM / (M-1))
* My formula is closer to yours than the AAVSO’s, but not quite the same.
* In my formula, SNR decreases as G gets larger. This must surely be correct? Higher gain cameras have more noise. In your formula, and the AAVSO’s, SNR increases with G. This doesn’t seem to make any sense.
* I think this formula behaves correctly when images are stacked. When you stack X images, sigma is divided by sqrt(X), and G is effectively divided by X. Net result: SNR increases as sqrt(X), as you expect.
* I’ve never done any aperture photometry in my life, so the above may or may not be complete nonsense.
Cheers,
Dominic
Thanks Dominic. Well beyond the call of duty to be doing stats at 2am!
That all looks reasonable. M/(M-1) and (M+1)/M aren’t quite the same but they are very similar and in any case tend to unity for the real situation of largish M. For the photon noise SG makes sense when G is defined as ADU/photon.
One further practical question about estimating G from the image. I am using the magnitude zero point derived by fitting an ensemble of stars to a catalogue. This links the magnitude (and hence photon flux) to the ADU count. Is there an alternative/better way of doing this?
I’ll try this later and will see if I can start getting some consistent SNR estimates.
One minor tweak….
You say that you need the outer annulus mean value and standard deviation. But, in cluttered fields, does using the median value help protect you from outliers?
Mean has better properties than median so I use an iterative approach which calculates the mean and rejects outliers.
It looks like gain is normally defined upside down (i.e. e-/ADU rather than ADU/e-) which is odd but would explain the reciprocal.
The one they use in VPHOT, described on page 4 here. Essentially you calculate the mean and SD then reject pixels more than 3 sigma from the mean and repeat until no more pixels are rejected. In the past I’ve used the median to get the sky level then used pixels below the sky level to determine the SD.
Grant,
Using the median isn’t a good way to exclude outliers, especially if a significant number of pixels are outliers (stars in a busy field). If, let’s say, 10% of pixels contain stars, then those outliers at the top of your distribution effectively mean that your median gives you the 55% percentile of the sky background distribution, which might be rather a noticeable skew towards the brighter end of the distribution.
Of course, that would be far better than taking the mean of a distribution containing stars, because the latter really would give you a nonsensical result. But to do this well, you really want to entirely get rid of the stars from your sample of pixels.
As Nick says, the way to do that is via “sigma clipping”. Take the mean and standard deviation of your pixel brightnesses. Exclude all pixels which are more than X standard deviations away from the mean. Iterate a few times.
I’ve come across this problem in the context of doing background subtraction from the images from my meteor cameras. In that context, the median has an additional issue that it’s quantised. If you take the median of a large sample of integer values between 0 and 255, you can’t get the noise below 1 ADU unit. With the mean of large numbers of samples, my sky background model sometimes reaches a formal precision of about ~ 0.2 ADU unit.
I would assume Nick has enough bit depth in his data that quantisation is unlikely to be a worry, though.
Cheers,
Dominic
If you really want to get rid of stars, DAOPHOT does a damn good job in my experience.
The work flow is to identify all the objects, then a set of unsaturated and unblended stars. From the latter you build a model of the point spread function. Note that the PSF may be position-dependent. Fit one or more PSF models to each detected object, removing ones which are “good enough” in some sense as you go. Repeat until you no longer have any “good enough” fits.
Thanks for the reply. I had wondered how much better sigma clipping might be. Looks easy enough to implement.
I have wondered if the skewness of the sample pixels could be calculated and that used as a threshold for pixel removal – though I suspect that might be more computationally intensive than sigma clip.
I shall have a play on some data I’m working on.
Just spotted this:
The reason I am doing this is that Astrometrica gets this wrong. In fact it reports two different SNRs in two different places
If by two different places you mean two different positions within the same image, it may not be wrong. As I noted in my earlier response the PSF can be position-dependent. The sky brightness will almost certainly be position-dependent (try measuring LBVs in the Andromeda galaxy!). Anyway, a fainter star (lower signal) will have a lower SNR than a brighter one even if the background is the same.
I am not quite understanding what you mean, in other words. I am but a bear of very little brain.
Astrometrica reports two different SNRs for the same measurement. In astrometry reports it uses something that Herbert Raab calls “peak SNR”. This is the SNR in the single brightest pixel. In the photometry report it uses what we would think of as the normal SNR definition (i.e. the SNR in the photometric aperture). Both I think use a process similar to what we have been discussing on this thread but I only really care about the photometric version since that one (or actually the log of it) is included in ADES astrometry that gets sent to the MPC. I’m trying to find out how Astrometrica estimates the gain, G (e-/ADU) from the image since the photometric SNR it reports does not scale correctly.
All the discussion about estimating the sky standard deviation is very interesting but for bright sources the photon noise dominates and getting G wrong really messes things up.
- You must be logged in to reply to this topic.