I. Introduction
In an excellent series of articles in CCD Astronomy (see note1), Michael Newberry discusses why and how we repair electronic artifacts in our digital CCD images. In particular, in MB5, he talks about flat fields and the "multiplicative" corrections they provide. He mentions using twilight sky flats, dome flats, and a hybrid of the two (illumination corrections). This webpage will dig a little farther into the two methods we use to flatten our images taken with the BGSU 0.5-m telescope and Ap6e CCD camera: twilight sky flats and illumination corrections.
In MN5, Newberry develops the idea that several causes alter light passing through a telescope in a multiplicative way (see note2). This includes vignetting, subtle dimming and shadowing within the optical system ("V"); dust on optical components leaving faint shadows ("s"); and variations in sensitivity to light from one pixel to the next due to variations in the quantum efficiency of the CCD material ("q") -- also see note3. He develops an equation that relates the signal S on the CCD (photons collected, or electrons measured) with the light falling on the primary mirror R and exposure time t:
S = R t V s q
We can visualize this equation by thinking about the different terms as functions of position across the CCD image, that is, S(x), R(x), V(x), s(x) and q(x). For example, we might consider the middle line of our CCD, Y=512, running across the image from X=1 to X=1024. The accompanying plot (click here for a permanent view of the image) shows some hypothetical functions that don't represent our actual CCD, but provide a clear and simple example.
R(x): In the lower left, we see the pattern of a perfectly uniform illumination of the primary mirror.
If the the telescope, camera and CCD were ideal and introduced no artifacts into the light, each pixel on the CCD would record 10,000 electrons. In general, this is not the case. The top panel of the figure shows four "response functions," or signatures, that describe the artifacts left in the image due to different imperfections in the optical system. They are:
V(x): The vignetting function varies slowly (has large-scale structure). The dip in the center is a very out of focus, partial shadow of the secondary mirror, and the dips at the edges may be due to partial shadowing by the filter wheel. Because the shadowing involves loss of light, the curve is everywhere < 1.0. A value of 0 would indicate no light reaches the CCD, a value of 0.95 means 5% of the light is lost, and a value > 1.0 would mean the light is somehow brighter than it might have been.
s(x): The dust shadowing function is also everywhere < 1.0, but it varies quickly in X (has small-scale structure) as we cut across a doughnut-shaped shadow "under" a dust speck located higher up the light path, on the filter or cameral window.
q(x): The average sensitivity value is often referred to as the CCD's quantum efficiency, the fraction of incoming photons that get converted into photoelectrons, typically 50-80% for modern CCDs. The sensitivity of the individual CCD pixels varies from one pixel to the next by a few percent around the average, producing the "jitter" in the figure.
II. Skyflats
When we take images of the clear twilight sky, we get a signal from the CCD that is the product of these functions:
Ss(x) = Rs(x) t V(x) s(x) q(x).
The plot in the bottom panel of our figure labeled "SkyFlat" is the product of the four functions in our example. Later in the night, we might take an image of a star cluster or other object that has a similar signal
So(x) = Ro(x) t V(x) s(x) q(x),
where Ro(x) is the function we wish to recover: the pattern of light incident on the telescope. We can do this by combining the observed So(x) with V(x), s(x), and q(x) from the sky flat and by assuming that R(x) is indeed uniform, or "flat" (this is why it is so important to take twilight skyflats when it is only perfectly clear and there is no significant moonlight; see note4). In practice, we do this by dividing the object image, pixel by pixel, by the flat field image:
So(x) / Ss(x) = [Ro(x) t V(x) s(x) q(x)] / [Rs(x) t V(x) s(x) q(x)] = Ro(x) / Rs(x)
The unknown response functions V(x), s(x) and q(x) cancel out (see note5). Since Rs(x) is constant or flat, we can replace it with simply Rs. Rearranging, we find Ro(x) = RsSo(x) / S(x)s. However, we don't know the value of Rs, a priori. If we set Rs to be the average value of S(x) across X (or the median, or mode), then Rs / Ss(x) is about 1.0, and the "flux" in So(x) is not greatly changed across the image. This is nice because it makes computing the signal-to-noise (see MN1 and MN2) and calculating the uncertainty in stellar magnitudes simpler.
The twilight skyflat method is the simplest and most direct method for flattening images, and is generally preferred if the sky is perfectly clear (no clouds or bright moon to give a non-uniform illumination) and if the observer can gather many skyflats, preferably during both dawn and dusk twilight.
III. Dome Flats
In many cases, it is not possible to get enough twilight skyflats in all the filters you are using to properly determine Ss(x), and hence use the twilight skyflat method. In MN5, Newberry describes how modern observatories have reflective screens, somewhat like a movie projector screen, mounted in the dome. During the afternoon, it is illuminated with color-balanced lamps (to simulate the spectrum of the night sky) and long sequences of "dome flats" with high signal-to-noise can be obtained at the observer's leisure.
However, these flats are rarely perfect either. Even with the best screens and lamps, it is difficult in practice to get a truly uniform illumination -- often there are variations at the 3-10% level. If we want to do stellar photometry with 1% accuracy, these flats will not be sufficient by themselves.
In our plot, the response function d(x) in the upper right describes how the illumination of the dome screen might deviate from a flat, uniform illumination -- the left edge receives about 8% less light than the average, and the right edge about 8% more, with a linear gradient in illumination across the CCD field. The plot in the bottom panel of our figure labeled "DomeF" is the product of the five functions in our example
Sd(x) = Rd(x) t V(x) s(x) q(x) d(x),
where Rd(x) is the idealized, uniform illumination we had hoped to shine on the telescope, and d(x) represents the amount by which the actual dome illumination deviates from this idealization (hence Rd(x) d(x) is the actual illumination received from the screen in the dome). There is a subtle but clear difference between the SkyFlat and DomeF functions in the plot owing to the d(x) term. Since this term is not known a priori and is difficult to measure, it is risky to use dome flats alone to flatten your data.
IV. Illumination Corrections
A more complicated process has been developed that combines the information from high signal-to-noise dome flats with low signal-to-noise twilight sky flats (or even multiple images of the dark sky during observations). Phil Massey (1997) provides a technical discussion of how this is done. In our implementation, we gather good skyflats from several nights over the course of the observing season -- images taken on very clear, moonless nights. We average them together to improve the signal-to-noise, reject stars and cosmic rays, and average over dusk and dawn flats that may have slightly different sky illumination. This master skyflat gives us the signal
Ss(x) = Rs(x) t V(x) s(x) q(x).
We also combine the domeflats taken on the night in question,
Sd(x) = Rd(x) t V(x) s(x) q(x) d(x).
Now we divide the skyflat by the domeflat to allow as many terms as possible to cancel,
Ss(x) / Sd(x) = [ Rs(x) t V(x) s(x) q(x) ] / [ Rd(x) t V(x) s(x) q(x) d(x) ] = Rs(x) / [Rd(x) d(x) ] = Rs / [Rd d(x) ]
where we have again replaced the R(x) functions with a constant equal to the average (median, or mode) across the chip. If we solve the equation for d(x), we have determined how much the dome illumination deviates from our ideal, flat pattern. However, the dust patterns s(x) may have changed somewhat between the nights the skyflats and domeflats were taken, so their cancellation may be imperfect. We therefore heavily smooth the d(x) function to remove these small-scale features and leave only the large-scale structure, d'(x). We call this the "illumination correction" function.
We can now apply this information to images of star clusters or other objects taken during the night that have a signal like
So(x) = Ro(x) t V(x) s(x) q(x).
We do this by dividing the observed image So(x) by the dome flat and then multiplying by the illumination correction d'(x):
So(x) / Sd(x) * d'(x) = [Ro(x) t V(x) s(x) q(x)] / [Rd(x) t V(x) s(x) q(x) d(x) ] * d'(x) = Ro(x) / Rd(x)
This time the dust patterns s(x) should be the same because the dome and object images were taken on the same night, so they should cancel more completely. The hypothetical dome response function d(x) should also cancel well with our measured illumination correction d'(x), provided our skyflats don't vary with time. Again assuming the hypothetical illumination Rd(x) is truly flat, we can replace it with a constant: the average (or median, or mode) across the chip, Rd. We then solve the equation for the function we wished to recover, the pattern of light from the object that was incident on the telescope,
Ro(x) = Rd d'(x) So(x) / Sd(x).
V. Summary
We have seen, from both graphical and mathematical perspectives, how sky flats can be used to remove the various artifacts left in a raw digital image. We extended the arguments to show how sky and dome flats can be combined to achieve the same goal for nights on which we were unable to obtain skyflats. In essence, the skyflats taken from other pristine observing nights, provide information on the large scale responses like vignetting, V(x), and dome illumination d(x), while the domeflats provide the small-scale structure due to dust rings, s(x), and pixel-to-pixel sensitivity variations, q(x). Since the dust features are most likely to change over times of hours or days as dust moves around on the optical elements and new dust settles out, it is important to take domeflats within a few hours of the object images you wish to flatten.
Before concluding, it is worth considering a few points. First, in our examples we have considered only one line of the image, leaving us relatively simple 1D expressions like R(x). In reality, IRAF handles all the lines of the image sequentially, giving us a 2D map of each function, like R(x,y). Second, each of the response functions may depend on the wavelength of the incident light, R(x). Therefore, we must take flat field images through each of the filters we plan to image with, and create separate master flat field images for each filter, so that we can properly flatten the images. Finally, because the skyflats were heavily smoothed in the process of building the illumination flat, we can employ much lower signal-to-noise images, and can often use fewer images that provides less complete rejection of stars, without badly degrading the signal-to-noise of our flattened object images (see MN2).
In reality, we can't achieve a perfectly flat distribution of incident light as shown as R(x) in our plot, because light streams in as individual photons, each with a random arrival time. As discussed in MN1, MN2 and our webage on signal to noise, this randomness is well described as noise that goes as the square root of the signal, N = sqrt(S). The function R'(x), shown at the lower right on our plot, contains this noise superimposed on our 10,000 electron signal from R(x). Each time we take an image, we get a different pattern of random noise, so in practice we take many images (flats) and average their results together, using a pixel rejection algorithm to detect and reject cosmic rays. This gives us a very high signal-to-noise master flat in which the noise is much smaller than the variations in any of the response functions V(x), s(x), q(x), and d(x) that we want to map and remove.
Endnotes
Note1: The five articles in CCD Astronomy by Michael V. Newberry, are The Signal to Noise Connection, 1994 Summer, p.34 (MN1); The Signal to Noise Connection Part II, 1994 Fall, p.12 (MN2); Recovering the Signal, 1995 Spring, p.18 (MN3), Dark Frames, 1995 Summer, p.12 (MN4); and Pursuing the Ideal Flat Field, 1996 Summer, p.18 (MN5).
Note2: The "additive" corrections, subtraction of the bias level image and subtraction of the dark current, have already been performed on the images we discuss here. The ordering of the steps is nicely described in MN3.
Note3: He also considers the effect of the camera shutter, which in the case of iris shutters like ours, leaves the center of the CCD exposed to slightly more light than the corners. We avoid this problem by avoiding exposure times shorter than 10 sec. Tests on our system suggest that the center of the chip is exposed to light for 0.1 sec (or less) longer then the corners. Thus, in a 10 sec exposure, the ratio is (tcenter - tcorner) / tcorner = (10.1 - 10.0 sec) / 10.0 sec = 0.01 or 1%. For longer exposures, it is less than 1%. Since 1% is the practical limit for obtaining good photometry due to lots of other factors, we set 10 sec as a safe lower-limit for exposure times at our telescope.
Note4: In theory, scattered moonlight should appear much like the scattered sunlight in the twilight sky, varying very slightly if at all across our small field of view. However, moonlight shining directly into the tube of our telescope (i.e., when the telescope is pointing within 90 deg of the moon) tends to reflect around and scatter onto the CCD in non-uniform ways. We have learned to be conservative and avoid taking skyflats when the moon is above the horizon.
Note5: In practice, we likely have different exposure times for the skyflat and object images, so they don't cancel out directly. This leaves a scaling factor that gets removed when we perform a photometric calibration on the instrumental magnitudes derived from the flattened images.