How Do I Perform CCD Image Reduction with Mathematica 4.2?

The purpose of this example is to demonstrate how Mathematica 4.2 can benefit from the support for data stored in the FITS format. FITS data often contain astronomical images from which can be extracted quantitative data. The data used in this example were obtained on May 9, 2000, in Flagstaff, Arizona, using the Lowell Observatory 31-inch telescope. The optical system included a liquid nitrogen-cooled CCD camera with an image resolution of 512 x 512 pixels.

All important system properties are stored in the header of each image. One exception is that the Filter property incorrectly reads "I" for infrared. In actuality, the data used here were taken using a V-band filter that has an effective wavelength of about 5,500 angstroms and a full width at half-max bandpass of about 890 angstroms.

Loading Reduction Frames

First we need to load the images necessary for the reduction process. The frames consist of a bias or zero frame, a flat-field frame, and the raw-image frame. In this case, the image was taken through a V-band filter, so the only flat-field frame needed is the one taken through the V-band filter. In all cases, ConversionOptions specifies that the header data should also be read in.

The bias and flat-field frames included in this section were created in another session from a series of frames. The average bias frame is the average of 10 separate bias frames. The average flat-field frame is an average of two separate flat-field frames taken in the V-band filter.

The bias frame is a measure of the "zero" point of the CCD chip, typically obtained by taking a zero-length exposure and simply reading out the pixel values from the chip. The bias frame is necessary so that when pixel counts are performed on the reduced image, they do not include the inherent electronic signal produced by the CCD itself. Since the variability of the bias signal is small across the CCD chip, a plot of the average bias would simply reveal a nearly uniform field that looks black. For this reason, the average bias frame will not be plotted here.

avgbias = Import["avgbias.fit",
    ConversionOptions -> {"Verbose" -> True}];

Next we read in the average flat-field frame and display it. The normalized flat-field frame is used to compensate for irregularities in the CCD sensitivity across the chip. These irregularities may be inherent to the chip itself or due to dust on the surface of the chip or on the filter being used. By exposing the CCD to a uniformly illuminated surface, the resulting image reveals these irregularities. By normalizing this field and then dividing the bias-corrected image by the normalized flat-field frame, the image can be corrected. Pixels that are less sensitive are corrected so that they are consistent with pixel values from more sensitive areas of the chip. Similarly, the more sensitive areas are corrected to be consistent with areas that are less sensitive.

avgvflat = Import["avgvflat.fit",
    ConversionOptions -> {"Verbose" -> True}];

ListDensityPlot[avgvflat[[1,2]],
    Mesh -> False, ImageSize -> 350];

The dark frame is a measure of the thermal noise produced by thermally excited electrons. In this case, the CCD chip was cooled with liquid nitrogen, so there was no appreciable thermal count. The dark frame was therefore omitted.

The image frame is an image of the giant elliptical galaxy M87. The raw image contains an overscan region that is 50 pixels wide on the right side. This region will be removed later. Plotting of the image frame is discussed in the next section.

image = Import["May08_0132.fit",
    ConversionOptions -> {"Verbose" -> True}];

We can also view the header, which contains important information such as the image type and the filter that was used.

TableForm[{"OBJECT", "BITPIX", "EXPTIME",
    "DATE-OBS", "OBSERVAT"}/.image[[1, 1]]]

Value -> 'M87'
Value -> 16.
Value -> 300.
Value -> '09/05/00'
Value -> 'lowell'

Reducing the Image

Reduction of an image is a process that processes a raw-image frame so that it is ready for useful analysis. In this case, the process involves subtracting a bias frame from all images and dividing the corrected image frame by a normalized flat-field frame.

We need to preserve the original header so that we can modify it to reflect the fact that pixels will be removed when the overscan region is removed.

oldheader = image[[1, 1]];

Since we will be removing the overscan region, we need to modify the header to indicate the reduction of the number of pixels.

newheader = oldheader/.("NAXIS1" -> _) ->
    ("NAXIS1" - > {"Value" -> 512., "Comment" -> ""});

We need to remove the overscan region from the flat-field frame and subtract the bias frame from it as well.

flatfielddata = Transpose[Take[Transpose[
    avgvflat[[1, 2]]], 512]] - avgbias[[1, 2]];

Similarly, we need to remove the overscan region from the raw-image frame.

olddata = Transpose[Take[Transpose[
    image[[1, 2]]], 512]];

To normalize the flat-field frame, we simply divide every pixel by the average pixel value of the entire frame.

normalizedflatfield = flatfielddata/
    ((Plus @@ Flatten @ flatfielddata)/(512*512.));

Now, we actually perform the reduction process on the image by subtracting the bias frame and then dividing by the normalized flat-field frame.

newdata = (olddata -
    avgbias[[1, 2]])/normalizedflatfield;

We finish by reconstructing the data into a header section and a data section. This is useful later if we wish to export the image to a new FITS image.

   workingimage = {{newheader, newdata}};

Plotting the Results

Now that the calculations have all been completed, we can see the results. First, we can view the raw-image frame. As the image shows, there are numerous "doughnuts" in the image plane that represent unfocused dust grains sitting in the light path. Also, the overscan region is visible on the right side of the image. Since the overscan region contains no image data, we removed it in the previous section.

ListDensityPlot[image[[1, 2]],
    Mesh -> False, ImageSize -> 350,
    AspectRatio -> Automatic];

Now, we see the reduced-image frame. The dust grains have been removed by the flat-fielding process, and the overscan region has been removed.

ListDensityPlot[workingimage[[1, 2]],
    Mesh -> False, ImageSize -> 350];

As a final step, we can adjust PlotRange to, in essence, adjust the contrast of the image. In this case, we can see that the core of M87 has a jet of material arising from a supermassive black hole at its core.

ListDensityPlot[workingimage[[1, 2]],
    PlotRange -> {75, 2500}, Mesh -> False, ImageSize -> 350];

The three FITS images used in this example are available for download.