I. Introduction:
In processing a set of CCD images, we often combine multiple images to improve the quality of the resulting image, specifically the ratio between the signal and the noise (or "signal-to-noise ratio," S/N). The signal represents information about the things in the image we care about -- stars, galaxies, nebulosity -- in the form of the number of counts (intensity) in the digital image at each pixel. The noise involves things such as the random noise incurred by reading out the CCD (read noise), and the random arrival of photons at the CCD (this noise source goes as the square root of the number photons received).
A second reason for combining multiple images is to detect and reject discrepant pixels. These may be pixels struck by cosmic rays, which arrive at random times and places on the CCD. Another source of discrepant pixels is ones that don't function properly, returning either no signal or a signal that is not proportional to the number of photons that fell on them ("bad pixels"). Other pixels may have high dark currents -- they produce a large number of thermal electrons per second unrelated to to the flux of photons falling on them -- that can not be corrected properly by our usual dark correction procedure ("hot pixels"). Fortunately, the number of such discrepant pixels tends to be small, though the number of cosmic ray hits and the severity of hot pixels do increase with the exposure time.
We use these techniques while combining multiple images of the sky -- a globular cluster or galaxy -- and when processing these images using the bias, dark, and flat-field images we took "off-sky." The accompanying plot (click here for a permanent view of the image) is specifically intended to represent the combining of bias frames taken with the BGSU Ap6e CCD, but the concepts are equally applicable to combining dark, flat, and "on-sky" object images, with our or any other CCD. For simplicity, let's consider only a small section of each image, a line of 99 pixels -- this might represent Y=429 and X=151 to 250 pix.
II. Individual Images
In panel (a) we see a plot of the intensity (in electrons, see note1) as a function of X-position along this line of pixels. The black line shows the counts from one image, the red line from a second image taken immediately after, and the blue line from a third image -- in each image we are considering the same line of pixels. The jagged up-down behavior centered on 4250 electrons intensity is due to the read noise of the CCD (see note2). The electronics can not read perfectly the number of electrons in each pixel, there is some measurement error, and since the reading from one pixel does not depend on the readings of surrounding pixels (we say the errors are uncorrelated), we get random scatter around the mean (see note3). The standard deviation of the points around the mean is a fair characterization of the size of this noise: for our CCD, it is about 12.7 electrons.
Notice that there are five cosmic rays affecting the data shown in panel (a). The sharp upward spikes at X=181, 238 and 246 in the black line, and at X=151 and 233 in the blue line, have high intensities because the single, high-energy cosmic ray that hit each of these pixels released many electrons. Since these spikes are caused by external events, we don't want to include their effects in our combined image -- we want to ignore, or "reject" them.
III. Averaging 3 Images and More
In order to reject the bad pixels and reduce the level of the noise, we can combine the images. At a given pixel, say X=190, we take the three intensity values from panel (a) -- for the black, red, and blue lines -- and average them: (4224 + 4254 + 4262)/3 = 4247. Doing this for every X in the line of pixels gives the curve shown in panel (b). Clearly, by averaging the uncorrelated errors in the individual images, we get a smoother curve. Colloquially, we "beat down the noise" by averaging multiple images.
Panel (c) shows the result of combining 15 different individual images, while panel (d) shows an image derived from 75 separate images. The more images we combine, the smaller the noise in the resulting image. This is why we strive to take ten or more bias and flat images every night at the telescope. We could do more, but the observing time increases linearly with the number of exposures (n) while the improvement to the combined image goes as sqrt(n), so there is a diminishing return. I think ten is a good compromise, while some people prefer 20 or 30. Notice that in panel (d), the noise is reduced enough so we can see a slight quadratic (concave-down) curvature in the intensity level as a function of X; this behavior is present in panels (a)-(c) but less obvious because it is masked by the noise. The bottom line is, by increasing the signal-to-noise ratio we improve our ability to see faint structures in our target.
IV. Rejecting Discrepant Pixels
This averaging process does not remove the effects of the cosmic rays. For example, at X=181, the average using three images would be (8452 + 4282 + 4240)/3 = 5658, well above the surrounding average. To detect and reject the effect of the cosmic ray, we employ a process called "sigma-clipping." We calculate the average of the intensities at each pixel and the standard deviation ("sigma"). At X=181, sigma = 1976. We then look at each of the intensities; if it is more than three sigma from the average (see note4), we reject the intensity and recompute the average and sigma based on the remaining values. In our example, the pixel hit by the cosmic ray is (8452-5658) = 2794 from the mean, or 1.4*sigma. Since it is not more than 3*sigma from the mean, we will not reject it. Notice that the spike still remains at this position in panel (b). If we had decided to reject the cosmic ray value, we would recompute the average and sigma from the remaining intensities: (4282 + 4240)/2 = 4261 and sigma=21. We could repeat this rejection procedure as many times as we liked (do 3 or 5 or 10 "iterations"), though this is only necessary if you are combining a large number of images and may have multiple cosmic rays falling on a given pixel.
This rejection method works better as the number of images gets larger, so we have more "good" values to draw on to define our average and standard deviation, thus making the cosmic ray a more extreme outlier. In panel (c), we have 15 images, and the initial mean and standard deviation are 4540 and 1082. The pixel hit by the cosmic ray is (8452-4540) = 3912 from the mean, or 3.6*sigma. This trips the sigma rejection, so we recompute the mean and standard deviation without the bright intensity value, giving 4260 and 10, respectively. All the cosmic rays that appeared in the individual images in panel (a) were successfully rejected in panel (c) when 15 images were used, and also in panel (d), when 75 images were used. There are other rejection methods available in IRAF (e.g., minmax, ccdclip, see note5) that may be successful when n is small.
V. Summary and Try-It-Yourself
So, as the number of images increases, both the signal-to-noise improves and the rejection of cosmic rays is more reliable, because our determination of the mean and standard deviation better reflect their true, underlying values (we get "better statistics"). The method above does not reject bad or hot pixels, however. If the pixel at X=213 were to be bad, it would be bad in every image, and the resulting mean and standard deviations would not help us detect and reject them. We may not be able to repair these pixels in the calibration images. However, when we are combining "on-sky" object images, say, of a globular cluster, we move the telescope a few arcsec between exposures so that the stars move across any bad or hot pixels (we "dither" the images). When we combine them, we shift the images so the stars are at a fixed pixel, and the hot/badpix are shifted around. The sigma-clipping method is then effective at detecting and rejecting hot and bad pixels as well as cosmic rays.
Here is an excel file with intensities from 20 different images, but at one pixel (perhaps X=184, Y=429). Use excel to recreate the sigma-clipping method described above to get a robust average intensity for this pixel (see note6). Iterate the process until no more pixels are rejected and your solution has "converged." When we run the combining tasks in IRAF, the computer does these calculations for each of the 1024*1024 pixels, completing each image in about a second -- ya gotta love computers!
More details about noise sources, how to calculate signal-to-noise ratios, and how the image processing steps impact signal-to-noise can be found in a nice set of articles by Michael Newberry in CCD Astronomy (1994, Summer, p.34; 1994 Fall, p.12; 1995 Spring, p.18; 1995 Summer, p.12; and 1996 Winter, p.18).
Finally, it is worth noting that we use this sigma-clipping algorythm in other places. For example, it is a nice way of removing the effect of discrepant data points when fitting a line or polynomial to a set of observed points: the CURFIT task in IRAF is one such place it can be used.
Endnotes
Note1: You may be familiar with the bias level on the BGSU CCD being around 500 adu. In the example here, I've converted that to electrons by multiplying by the known gain of the CCD: (500 adu) * (8.5 electrons/adu) = 4250 electrons. Since one photon frees one electron and the Poisson noise is proportional to the square root of the number of photons falling on a pixel, it is appropriate to do the noise calculations in units of electrons, not adu.
Note2: In each of the panels (a)-(d), the dotted line with an intensity of 4250 electrons is drawn to guide the eye. The y-axis scale of each of the panels is the same; the numbers were removed from panels (b) and (c) for clarity.
Note3: If this represented an "on-sky" object image, the jagged behavior would include a contribution from the background sky brightness (whether due to light pollution, moonlight, natural skyglow, or background objects like nebulosity or unresolved galaxies or stars) in addition to the read noise.
Note4: The choice to reject intensities more than 3*sigma from the mean is a bit arbitrary. We could choose 5*sigma if we wanted to be more conservative, or 2*sigma. If we assume the noise has a normal (Gaussian) distribution about the mean, we will accidentally reject the rare "noise spikes" about 0.3% of the time when doing 3*sigma rejection, but would expect to reject about 5% of these spikes when doing 2*sigma rejection. They would be incredibly rare in 5*sigma rejection, hence it is more conservative.
Note5: The option "minmax" removes the brightest n values from the stack for a given pixel, along with the faintest m values (the user selects n & m) and averages the remaining values; its a quick but imperfect way of removing cosmic rays. The option "ccdclip" predicts what sigma should be based on the known readnoise of the CCD, the signal (N=sqrt(S)), and other factors, rather than the observed sigma from the pixel values themselves. A cosmic ray hit would contribute to the latter, but not the former, so the cosmic ray value would tend to be many sigma from the mean and easily detected/rejected.
Note6: The excel functions average(z1:z2), stdev(z1:z2) and abs(z1:z2) may help you.