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