Dark current: the ideal case
Contents
3.1. Dark current: the ideal case#
import numpy as np
from scipy import stats
%matplotlib inline
from matplotlib import pyplot as plt
from image_sim import dark_current, read_noise
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[1], line 7
4 get_ipython().run_line_magic('matplotlib', 'inline')
5 from matplotlib import pyplot as plt
----> 7 from image_sim import dark_current, read_noise
File ~/work/ccd-reduction-and-photometry-guide/ccd-reduction-and-photometry-guide/notebooks/image_sim.py:6
3 import numpy as np
5 from astropy.modeling.models import Gaussian2D, RickerWavelet2D, Const2D
----> 6 from photutils.datasets import (make_random_gaussians_table,
7 make_gaussian_sources_image)
8 from photutils.aperture import EllipticalAperture
10 # To use a seed, set it in the environment. Useful for minimizing changes when
11 # publishing the book.
ImportError: cannot import name 'make_random_gaussians_table' from 'photutils.datasets' (/usr/share/miniconda/envs/test/lib/python3.11/site-packages/photutils/datasets/__init__.py)
# Use custom style for larger fonts and figures
plt.style.use('guide.mplstyle')
3.1.1. A dark frame measures dark current#
Recall that dark current refers to counts (electrons) generated in a pixel because an electron in the pixel happens to have enough energy to “break free” and register as a count. The distribution of electron thermal energies in a pixel follows a Maxwell-Boltzmann distribution in which most electrons have energy around \(kT\), where \(T\) is the temperature of the sensor and \(k\) is the Boltzmann constant. There is a distribution of energies, though, and occasionally an electron will be high energy enough to jump to the conducting band in the chip, registering the same as an electron excited by a photon. Since the Maxwell-Boltzmann distribution depends on temperature, the rate at which dark current appears in a pixel is also expected to depend on temperature.
A dark frame (also called a dark image) is an image taken with your camera with the shutter closed. It is the sum of the bias level of your camera, the readout noise, and the dark current.
You measure the dark current in your camera by taking dark frames.
3.1.1.1. Convenience functions#
A couple of functions that will be reused a few times in this notebook are defined below.
A function to generate dark frames
This function generates a set of dark frames. It will be used several times below.
def generate_dark_frames(n_frames,
dark_rate=0,
noise=0,
gain=1,
exposure=10,
image_size=500):
base_image = np.zeros([image_size, image_size])
darks = np.zeros([n_frames, image_size, image_size])
for n in range(n_images):
darks[n] = dark_current(base_image, dark_rate, exposure, gain=gain, hot_pixels=False)
if noise > 0:
darks[n] = darks[n] + read_noise(base_image, noise, gain=gain)
return darks
A function to plot the distribution of dark counts
The function below plots the:
def plot_dark_with_distributions(image, rn, dark_rate,
n_images=1,
exposure=1,
gain=1,
show_poisson=True,
show_gaussian=True):
"""
Plot the distribution of dark pixel values, optionally overplotting the expected Poisson and
normal distributions corresponding to dark current only or read noise only.
Parameters
----------
image : numpy array
Dark frame to histogram.
rn : float
The read noise, in electrons.
dark_rate : float
The dark current in electrons/sec/pixel.
n_images : float, optional
If the image is formed from the average of some number of dark frames then
the resulting Poisson distribution depends on the number of images, as does the
expected standard deviation of the Gaussian.
exposure : float
Exposure time, in seconds.
gain : float, optional
Gain of the camera, in electron/ADU.
show_poisson : bool, optional
If ``True``, overplot a Poisson distribution with mean equal to the expected dark
counts for the number of images.
show_gaussian : bool, optional
If ``True``, overplot a normal distribution with mean equal to the expected dark
counts and standard deviation equal to the read noise, scaled as appropiate for
the number of images.
"""
h = plt.hist(image.flatten(), bins=20, align='mid',
density=True, label="Dark frame");
bins = h[1]
expected_mean_dark = dark_rate * exposure / gain
pois = stats.poisson(expected_mean_dark * n_images)
pois_x = np.arange(0, 300, 1)
new_area = np.sum(1/n_images * pois.pmf(pois_x))
if show_poisson:
plt.plot(pois_x / n_images, pois.pmf(pois_x) / new_area,
label="Poisson dsitribution, mean of {:5.2f} counts".format(expected_mean_dark))
if show_gaussian:
# The expected width of the Gaussian depends on the number of images.
expected_scale = rn / gain * np.sqrt(n_images)
# Mean value is same as for the Poisson distribution
expected_mean = expected_mean_dark * n_images
gauss = stats.norm(loc=expected_mean, scale=expected_scale)
gauss_x = np.linspace(expected_mean - 5 * expected_scale,
expected_mean + 5 * expected_scale,
num=100)
plt.plot(gauss_x / n_images, gauss.pdf(gauss_x) * n_images, label='Gaussian, standard dev is read noise in counts')
plt.xlabel("Dark counts in {} sec exposure".format(exposure))
plt.ylabel("Fraction of pixels (area normalized to 1)")
plt.grid()
plt.legend()
3.1.2. Dark current theory#
The expected signal in a dark frame exposure of time \(t\) is proportional to \(t\). If we call the dark electrons in an exposure \(d_e(t)\) and the dark current \(d_c(T)\), where \(T\) is the temperature, then
For liquid-cooled cameras, particularly ones cooled by liquid nitrogen, the operating temperature doesn’t change. For thermoelectrically-cooled cameras you are able to set the desired operating temperature. As a result, you should be able to ignore the temperature dependence of the dark current.
The thermoelectric coolers can usually cool by some fixed amount below the ambient temperature. Though in principle you could choose to always cool by the same fixed amount, like \(50^\circ\)C below the ambient temperature, there is an advantage to always running your camera at the same temperature: dark frames taken on one date are potentially useful on another date. If the operating temperature varies then you need to make sure to take dark frames every time you observe unless you carefully characterize the temperature dependence of your dark current.
It will turn out that for practical reasons — not all pixels in your camera have the same dark current — it is usually best to take dark frames every time you observe anyway.
3.1.3. Illustration with dark current only#
For the purposes of illustrating some of the properties of dark current and dark frames, we’ll generate some simulated images in which the counts are due to dark current alone. We’ll use these values:
Dark current is \(d_c(T) = 0.1 e^-\)/pixel/sec
Gain is \(g = 1.5 e^-\)/ADU
Read noise is 0 \(e^-\)
dark_rate = 0.1
gain = 1.5
read_noise_electrons = 0
3.1.3.1. Dark current is a random process#
The dark counts in a dark frame are counts and so they follow a Poisson distribution. The plot below shows the dark current in a number of randomly chosen pixels in 20 different simulated images each with exposure time 100 sec. Note that the counts vary from image to image but that the average is very close to the expected value.
The expected value of the dark counts for this image is \(d_e(t)/g = 6.67~\)counts.
exposure = 100
n_images = 20
n_pixels = 10
image_size = 500
plt.figure(figsize=(20, 10))
darks = generate_dark_frames(n_images, dark_rate=dark_rate, noise=read_noise_electrons,
gain=gain, exposure=100, image_size=image_size)
# Pick several random pixels to plot
pixels = np.random.randint(50, high=190, size=n_pixels)
pixel_values = np.zeros(n_images)
pixel_averages = np.zeros(n_images)
for pixel in pixels:
for n in range(n_images):
pixel_values[n] = darks[n, pixel, pixel]
plt.plot(pixel_values, label='pixel [{0}, {0}]'.format(pixel), alpha=0.5)
pixel_averages += pixel_values
plt.plot(pixel_averages / n_pixels,
linewidth=3,
label='Average over {} pixels'.format(n_pixels))
# plt.xlim(0, n_images - 1)
plt.hlines(dark_rate * exposure / gain, *plt.xlim(),
linewidth=3,
label="Expected counts")
plt.xlabel('Image number')
plt.ylabel('Counts due to dark current')
plt.legend()
plt.grid()
3.1.3.2. Dark counts follow a Poisson distribution#
The distribution below shows a normalized histogram of number of pixels as a function of dark counts in each pixel for one of the simulated dark frames. Overlaid on the histogram is a Poisson distribution with a mean of \(d_e(t_{exp}) = d_C(T) * t_{exp} / g\), where \(t_{exp}\) is the exposure time.
plt.figure(figsize=(20, 10))
h = plt.hist(darks[-1].flatten(), bins=20, align='mid', density=True,
label="Histogram of dark current counts in single image");
bins = h[1]
pois = stats.poisson(dark_rate * exposure / gain)
pois_x = np.arange(0, 20, 1)
plt.plot(pois_x, pois.pmf(pois_x),
label="Poisson dsitribution, mean of {:5.2f} counts".format(dark_rate * exposure / gain))
plt.xlabel("Dark counts in {} exposure".format(exposure))
plt.ylabel("Number of pixels (area normalized to 1)")
plt.legend()
plt.grid()
3.1.4. Illustration with dark current and read noise#
Now let’s run through the same couple of plots with a non-zero read noise. For the sake of illustration, we’ll look at two cases:
Moderate read noise of 10 \(e^-\) per read, typical of a low-end research-grade CCD
Low read noise of 1 \(e^-\) per read
In both cases we’ll continue with the parameters above to generate our frames:
Dark current is \(d_c(T) = 0.1 e^-\)/pixel/sec
Gain is \(g = 1.5 e^-\)/ADU
Exposure time 100 sec
With those choices the expected dark count is 6.67 count, which is 10 \(e^-\). That is, not coincidentally, one of the values for read noise that was chosen.
The dark current rate in this example is somewhat higher than that of a real CCD. A large dark current is used here to demonstrate the interaction of dark current and read noise, and the exposure time has been chosen that expected dark electrons is 10\(e-\).
3.1.4.1. High noise#
In this first case, the read noise and the dark current are both 10\(e^-\). As we will see shortly, the dark frame really measures read noise, not dark current, under these circumstances.
high_read_noise = 10
darks = generate_dark_frames(n_images, dark_rate=dark_rate, noise=high_read_noise,
gain=gain, exposure=exposure, image_size=image_size)
image_average = darks.mean(axis=0)
plt.figure(figsize=(20, 10))
plot_dark_with_distributions(darks[-1], high_read_noise, dark_rate,
n_images=1, exposure=exposure, gain=gain)
plt.xlim(-20, 30);
This dark frame measures noise, not dark current
The pixel distribution is clearly a Gaussian distribution with a width determined by the read noise, not the underlying Poisson distribution that a dark frame is trying to measure. The only way around this (assuming the dark current is large enough that it needs to be subtracted at all) is to make the exposure long enough that the expected counts exceed the dark current.
We explore that case below by adding in a much smaller amount of noise.
Combining does not recover the dark current
The plots below show the result of averaging several (20 in this case) dark frames together. We might hope that averaging reduces the noise enough that the distribution of pixels becomes the Poisson distribution expected for a dark frame. It does not.
Note: To correctly calculate the expected distribution of counts in the average of the dark frames you should start with the expectations for the sum of the images, then scale counts down by the number of images and the distribution up by the same factor so that the distribution remains normalized.
plt.figure(figsize=(20, 10))
plot_dark_with_distributions(image_average, high_read_noise, dark_rate,
n_images=20, exposure=exposure, gain=gain)
plt.title("Distribution of counts for the average of 20 dark images, high noise case");
Even the combination of 20 dark frames doesn’t lead to a Poisson distribution; the distribution is still a Gaussian whose width is determined by the read noise and the number of images.
3.1.4.2. Low noise#
In this case the read noise is 1 \(e^-\), lower by a factor of ten than the expected dark current for this exposure time, 10\(e^-\).
This case almost never occurs in real CCDs, but it illustrates that with a low enough read noise you can recover the Poisson distribution of the dark current.
low_read_noise = 1
darks = generate_dark_frames(n_images, dark_rate=dark_rate, noise=low_read_noise,
gain=gain, exposure=exposure, image_size=image_size)
image_average = darks.mean(axis=0)
plt.figure(figsize=(20, 10))
plot_dark_with_distributions(darks[-1], low_read_noise, dark_rate,
n_images=1, exposure=exposure, gain=gain)
#plt.ylim(0, 0.8)
plt.xlim(-10, 20)
plt.title("Distribution of counts for the average of one dark images, low noise case");
plt.figure(figsize=(20, 10))
plot_dark_with_distributions(image_average, low_read_noise, dark_rate,
n_images=20, exposure=exposure, gain=gain)
plt.ylim(0, 0.8)
plt.xlim(0, 14)
plt.title("Distribution of counts for the average of 20 dark images, low noise case");
Dark frame is measuring dark current
Unlike the high noise case, the pixels in a single low-noise dark frame follow the Poisson distribution expected for the dark current in this simulated camera. When multiple images are combined the pixels still follow the Poisson distribution, as expected in this case.