ProFit: Sigma Maps
In ProFit: Fit Projected 2D Profiles to Galaxy Images

Get the latest version of ProFound and ProFit:

library(devtools)
install_github('asgr/ProFound')
install_github('ICRAR/ProFit')

Set global evaluate (basically TRUE for GitHub version and FALSE for CRAN):

evalglobal=FALSE

First load the libraries we need:

library(ProFit)
library(ProFound)

The authors of ProFit have answered many (many) questions regarding what sigma maps mean / represent, and how to compute them in practice for input to profitSetupData. This short vignette will clarify both of these questions.

Make some fake data:

modellist=list(
  sersic=list(
    xcen=c(200, 200),
    ycen=c(200, 200),
    mag=c(17, 13),
    re=c(5, 10),
    nser=c(4, 1),
    ang=c(0, 135),
    axrat=c(1, 0.3),
    box=c(0,0)
  )
)

psf=profitMakeGaussianPSF(fwhm=3, dim=c(25,25))

image_elec=profitMakeModel(modellist=modellist, psf=psf, dim=c(400, 400), magzero=30)$z

Here we use a magzero point of 30, since this is a pretty typical value for astronomical Sloan r-band images that are in counts when we are working in the AB magntiude system (it is beyond the scope of this document to tell you why, but in short it is because galaxies span a fairly narrow range of surface brightness and CCDs tend to have well depths <100k electrons).

Now we can check the image and the counts:

magimage(image_elec)
maghist(image_elec, log='x', majorn=c(10,5))

Notice that the image counts tail off at about 36k. Given that typical CCDs have well depths <100k and are most linear in the middle regime of electron counts, this image represents an ideally exposed image.

An important point to note is that the image is represented, in its rawest form, in photo-electron counts (not photon counts). For modern CCDs the difference tends not to be much (quantum efficiency / QE = 0.9 typically in the r-band). In any case, we care about the number of electrons excited by photons, not the number of photons that do the exciting. Exactly why is a more detailed explanation, but in brief QE is a stochastic binomial filtering process, not an exact fraction. Since every photon has, say, a 50% chance of either exciting an electron or not, the net result is that it is the photo-electrons that produce the image shot-noise and have a Poisson distribution (error goes as $\sqrt{N}$ for large $N$). A quick R proof:

RanPois=rpois(n=1e4, lambda=400)
RanPois_QE=rbinom(n=1e4, size=RanPois, prob=0.5)
RanPois_PureFilter=RanPois/2

sd(RanPois_QE)
sd(RanPois_PureFilter)

The proper QE process has $\sqrt{2}$ more noise, and is the correct type of noise to consider for a sigma map.

Since the image we have made is in photo-electron counts by construction, we should formally use Poisson statistics to represent the uncertainty per pixel. In practice in the regime that matters ($N>20$) we can approximate the behaviour with Normal statistics. So now we can make our object only shot noise sigma map:

sigma_elec=sqrt(image_elec)

To our image we should add a realistic component of sky flux/counts. This clearly introduces its own component of image shot-noise. Typically sky photon-electron counts tend to number 100s and sky noise tends to be in the regime of 10s of counts for well exposed data (e.g. $\sqrt{400}=20$). To add the error correctly we have to add things together in quadrature.

sky_elec=400
sky_elec_noise=sqrt(sky_elec)

image_elec=image_elec+sky_elec
sigma_elec=sqrt(sigma_elec^2+sky_elec_noise^2)

Real data also tends to have a small degree of dark current (a few counts typically), which brings its own shot, and also bias counts and the associated read-noise. The dark current follows Posisson count statistics, but the read-noise tends to be less than you might expect from the bias electron counts. This gets into complicated CCD issues, but basically you are adding in noise from various imperfect column reading and analogue to digital data conversion process, as well as some noise from the artifical addition of electrons to each pixel (these are not generated by true Poisson processes). Here we will assume 9 counts of dark current producing 3 counts noise, and 900 counts of bias that have 10 counts of read-noise (less than $\sqrt{900}=30$).

dark_elec=9
dark_elec_noise=sqrt(dark_elec)

bias_elec=900
read_elec_noise=10

image_elec=image_elec+dark_elec+bias_elec
sigma_elec=sqrt(sigma_elec^2+dark_elec_noise^2+read_elec_noise^2)

Already we can see that we cannot trivially go from image electron counts to the sigma electron map:

maghist(sqrt(image_elec)/sigma_elec)

Things get more complicated when dealing with real astronomical data, beacuse you tend to be provided with bias (always), dark (always) and sky (often) subtracted images. Let's say they do the job perfectly:

image_elec_sub=image_elec - sky_elec - dark_elec - bias_elec

Next astronomers tend to put images into different Astronomical Data Unit (ADU) schemes where they want one ADU to represent 0 magnitudes. In this case the magzero point input to ProFit is 30, and the image pixels become manipulated as such:

magzero=30
gain_elec_ADU=1/10^(-0.4*magzero)
image_ADU_sub=image_elec_sub/gain_elec_ADU

The gain is nearly always expressed in this sense (photo-electrons per ADU) rather than ADU per photo-electron, but I have seen both used informally (ouch). Just to be clear (since it is confusing) there are at least two types of gain in astronomy:

Analogue to digital conversion (ADC) gain applied at the CCD level when reading off (usually a number between 0.5-2). It is really there to get counts close to a power of 2 (since the readout will always be $0 \rightarrow 2^n-1$) given the well depth (which can be an awkward number like 100,000 electrons).
Gain that takes you from arbitrarily scaled ADUs to photo-electrons.

Both of these types of gain have the same formal units of photo-electrons per ADU (so photo-electrons = ADUs x gain). Whether or not the origin of the gain is from an on-chip multiplication process or a later scaling when changing the zero-point the aim is the same: it is a number that can take us all the way from arbitrary pixel units to the true number of excited electrons in each pixel well. This is vital, since it is these excited electrons that behave with true Poisson statistics.

To this exact image (remember, this is our exact model still) we can add the various noise terms back in using a Normal sampling process. For increased accuracy we could use a Poisson process, but since the per pixel counts are fairly large ($N>20$), this is not necessary.

image_ADU_sub_noise=image_ADU_sub+rnorm(length(image_ADU_sub), mean=0, sd=sigma_elec)/gain_elec_ADU

magimage(image_ADU_sub_noise)

Notice this fairly realistic image is much shallower than the first one we made. We cannot observe so far into the low surface brightness wings of the galaxy due to the various noise components that do not come from the galaxy object shot noise (read, dark, sky). This is why it is important to observe galaxies in dark sites- galaxies are pretty dim compared to the typical night sky we observe on the ground.

Transform data back to a sigma map

The question is, given the type of data in image_sub_ADU_noise and knowledge that the gain_e_ADU=1e12, how do you calculate the sigma image when this has not been explicitly provided? In pratice there are two routes (both of which are imperfect).

In one you are told the sky pedestal subtracted (in ADUs), the read noise (in photo-electrons), the dark noise (in photo-electrons) and the gain (in photo-electrons per ADU):

$$\frac{\sqrt{(image[ADU]+sky[ADU]).gain+\sigma_{read}[elec]^2+\sigma_{dark}[elec]^2}}{gain}$$

So in practice, with the data available, we would approximate this from image_sub_ADU_noise with:

sky_ADU=sky_elec/gain_elec_ADU
sigma_est_1=sqrt((image_ADU_sub_noise+sky_ADU)*gain_elec_ADU + read_elec_noise^2 + dark_elec_noise^2)/gain_elec_ADU

Now we can compare this estimated sigma map against the exact answer we know from above. It should be a delta function at 1 if it was precisely correct everywhere on the image.

maghist((sigma_est_1/sigma_elec)*gain_elec_ADU)

So even with this high quality of meta information we find a few percent uncertainty in our estimated sigma map. Even this scenario is a bit optimisitic in practice, but it is achievable with data sources that include a lot of meta data and helpful online documentation (e.g.\ SDSS, HST).

More realisticially, you will need to use the gain to estimate the shot noise of your object, and estimate the other noise terms from the noise present in the image itself.

First you will need a reasonable segmentation map (of the type you would pass to segim in profitSetupData). Here we will be very conservative and only use the faintest 50% of pixels from the original (noise free) image we made:

segim=image_elec>=quantile(image_elec,0.5)

This will look like a white circle when we plot it

magimage(segim)

Now we can look at the sky data in the non-object pixels:

  sky_pix_ADU=image_ADU_sub_noise[segim==0]
maghist(sky_pix_ADU)

You should check how reasonable our selection of sky pixels is by comparing the negative sky pixels to their positive counterparts (remember, we have an image with subtracted sky so the noise should be symmetric around 0). There are no negative flux astronomical sources, so you tend to see a slight excess for positive counts if the masking is not agressive enough. In general more agressive masking of sources is wise when trying to estimate noise in this manner.

magplot(density(log10(sky_pix_ADU[sky_pix_ADU>0]), bw=0.01), col='red', type='l')
lines(density(log10(abs(sky_pix_ADU[sky_pix_ADU<0])), bw=0.01), col='blue')

Luckily for us, our segim mask has worked pretty well, and we get nice symmetrical values either side of 0. This means we can easily estimate the noise. A nice side effect is that the standard-deviation of sky pixels contains the sky shot-noise, read-noise and dark-noise, i.e.:

$$\sigma_{skypix}[ADU]=\sqrt{\sigma_{sky}[ADU]^2+\sigma_{read}[ADU]^2+\sigma_{dark}[ADU]^2}$$

sky_pix_ADU_noise=sd(sky_pix_ADU)
sigma_est_2=sqrt(image_ADU_sub_noise/gain_elec_ADU+sky_pix_ADU_noise^2)

maghist(sigma_est_2/sigma_elec*gain_elec_ADU)

As a convenience there also exists a helper function profitMakeSigma that does the scaling and error addition for you (if you are feeling lazy, or want to sanity check your own thought process). Below we compare both of the above methods to estimate the sigma maps:

sigma_est_3=profoundMakeSigma(image=image_ADU_sub_noise, sky=0, skyRMS=sky_elec_noise, gain=gain_elec_ADU, readRMS=read_elec_noise, darkRMS=dark_elec_noise, sky_units='elec', read_units='elec', dark_units='elec')
magimage(sigma_est_3)

sigma_est_4=profoundMakeSigma(image=image_ADU_sub_noise, sky=0, skyRMS=sky_pix_ADU_noise, gain=gain_elec_ADU)
magimage(sigma_est_4)

The only difference is the helper function does not allow negative sky counts to reduce the skyRMS floor:

magplot(sigma_est_2, sigma_est_3, pch='.', log='xy', grid=TRUE)

The outputs of using profitMakeSigma either way is almost identical:

magplot(sigma_est_3, sigma_est_4, pch='.', log='xy', grid=TRUE)

We can also use the segmentation map we made before to identify pure sky areas where the error should simply be the skyRMS:

sigma_est_5=profoundMakeSigma(image=image_ADU_sub_noise, objects=segim>0, sky=0, skyRMS=sky_pix_ADU_noise, gain=gain_elec_ADU)
magimage(sigma_est_5)

magplot(sigma_est_2, sigma_est_5,pch='.',log='xy', grid=TRUE)
abline(h=sky_pix_ADU_noise,col='red')

The histogram a few plots above looks much like our first estimate, but it is a tiny bit skewed to lower values. In practice both produce decent Normal residuals (sd~1) when we look at (Data-Model)/Sigma. The best we can possibly do is:

sd((image_ADU_sub_noise-image_ADU_sub)/sigma_elec)*gain_elec_ADU

Let us see how our approximate sigma maps compare. Note if our sigma estimates are good then the sigma residual should look like featureless noise:

magimage((image_ADU_sub_noise-image_ADU_sub)/sigma_est_1)
magimage((image_ADU_sub_noise-image_ADU_sub)/sigma_est_2)

magplot(density((image_ADU_sub_noise-image_ADU_sub)/sigma_est_1))
curve(dnorm, col='red', add=T)
magplot(density((image_ADU_sub_noise-image_ADU_sub)/sigma_est_2))
curve(dnorm, col='red', add=T)

sd((image_ADU_sub_noise-image_ADU_sub)/sigma_est_1)
sd((image_ADU_sub_noise-image_ADU_sub)/sigma_est_2)

We can check how sensitive this approach might be to error:

magimage((image_ADU_sub_noise*1.01-image_ADU_sub)/sigma_est_1)

Basically, both our two estimates are no worse than the correct answer. It is worth considering why there is a difference. Ultimately the sigma map should be derived from the exact model, not a noisy sample from it. I.e. we might be observing a pixel that intrinsically (averaged over infinite time and images) produces, e.g., 10,000 / 100 / 10 photon-electron counts in the time we integrate our image for, but in any given image we do not observe the exact answer. How inaccurate is this assumption?

maghist(sqrt(rpois(1e3,1e4)), breaks=20)
maghist(sqrt(rpois(1e3,1e2)), breaks=20)
maghist(sqrt(rpois(1e3,1e1)), breaks=20)

This is not too bad in general (when photon-electron counts are high). It means if our source is this bright we should expect to be able to estimate the error for pixels with 10,000 / 100 / 10 photon-electron counts to better than 1% / 5% / 20% accuracy. In fact there is a funny effect that occurs with Poisson statistics that makes computing this very easy- taking the square-root of a Poisson distribution is variance stabilising, so for $N>5$ the square-root of the Poisson distribution gives a standard-deviation close to 0.5, thus the approximate relative error in the above assumption simply becomes $0.5/\sqrt{cnts_{PE}}$.

For details on this look at the article on Anscombe transforms on Wikipedia, but in brief the transform $A : x \rightarrow 2\sqrt{x + 3/8}$ transforms all Poisson distributions so they have a mean $\mu_{new}=2\sqrt{\mu_{old}+3/8}-1/(4\sqrt{\mu_{old}}))$ and standard-deviation $\sigma_{new}=1$. Hence our simpler tranform (without the additive $3/8$ and factor 2 multiplier) generates a standard-deviation uncertainty close to 0.5 for all pixels when estimating their true mean.

That rounds off this vignette. You should now have a good idea how to approch making a sigma map for feeding into profitSetupData.