7 Image correction
7.1 Introduction
The number of counts recorded in the image on a pixel is a linear function of the flux falling on that pixel in the CCD:
Both a and b vary by pixel, and may also vary be time. The time variation is such that it is generally wise to recalculate them on a daily basis, although, on the PT, in practice one can get away with using values from nearby nights.
The value of b, the bais, can be measured using a combination of zero second exposures and overscans. The linear constant, a, can the measured using observations of a uniformly illuminated screen, or perhaps a blank area of the sky.
The process described below is neigther the optimal or most convenient appreach. Rather, it is designed to clarify what needs to be done.
This page only describes the most common image correction operations: bias subtraction and flat fielding. A few others are sometimes needed. Charge Transfer Efficiency (CTE) correction accomodates inefficiency in the transfer of charge between pixels durning readout. Dark correction removes the natural accumulation of charge over time not due to hits by photons. The former is unimportant for most uses of PT data, and the camera on the PT has a negligible dark current, which can be ignored.
7.2 Bias subtraction
7.2.1 The Master Bias
The sample data directory contains a sequence of bias images, exposures 114031 through 114041. These are all zero second exposures, readouts of the CCD taken without opening the shutter. These images may still, however, contain cosmic rays. A quick-and-dirty way to get a bias frame is to median filter the images. I will use Python and a few Python libraries: pyfits, a Python FITS reader; numarray, which handles arrays such as FITS images; and numdisplay, which provides a Python interface to ds9 for veiwing data arrays.
Start python as in Python, then
Python 2.5.1 (r251:54863, Dec 20 2007, 12:03:19)
[GCC 3.2.3 20030502 (Red Hat Linux 3.2.3-59)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import pyfits
>>> from numpy import *
>>> import numdisplay
>>> biasData = array([pyfits.getdata("mdR-%08d.fits" % expId) for expId in range(114031,114042)])
>>> bias = median(biasData)
If we want, we can save our bias as a FITS file:
>>> hdu=pyfits.PrimaryHDU(bias)
>>> hdu.header.add_comment("Median of PT exposures 114031 through 114042")
>>> pyfits.writeto("bias.fits",bias)
We can display the FITS image using ds9. Alternately, the numdisplay package provides a Python interface to ds9 for veiwing data arrays. We do not even have to save the array to use this interface.
>>> median(bias) array([ 1346., 1346., 1346., ..., 1272., 1272., 1272.], dtype=float32) >>> numdisplay.display(bias,z1=1200,z2=1400) Image displayed with Z1: 1200 Z2: 1400The numdisplay display quantizes the image to 256 gray levels before sending it to ds9, and does not always choose useful transformations automatically. It has a number of parameters that let you control this. In the above example, I use the median function to get the medians of the columns in the bias array to get an idea of what the bias level, and force the min (z1) and max (z2) limits of the transformation accordingly.
Now, load the data image and apply the bias, and take a look:
>>> rawImage = pyfits.getdata("mdR-00114174.fits")
>>> dbImage = rawImage - bias
>>> numdisplay.display(dbImage)
Image displayed with Z1: -6.0 Z2: 64182.0
>>> numdisplay.display(dbImage,z1=0,z2=1000)
Image displayed with Z1: 0 Z2: 1000
Here, I chose the scaling using an educated guess.
7.2.2 The Overscan Region
While the master bias frame takes into account pixel to pixel differences in bais, the bias can change significantly with time over the course of the night; the true bias of the data image will be slightly different. These differences can be estimated using overscan regions: pixels in the final image that have been read through the CCD amplifier but that were not part of the exposed images. For the PT, these are data values are obtained from the CCD by reading more pixels per row and more rows per readout than are actually present on the CCD. In other instruments, actual obscured pixels are used.
There are no officially standard FITS keywords for designating what sections of the a readout are overscan, but the information is generally present somewhere in the header. The PT data uses the header keywords BIASSEC0, BIASSEC1, DATASEC0, and DATASEC1 to designate the areas of the CCD used for the row bias and data sections for amplifiers 1 and 2, respectively. There is also a column bias section, but it is not explicitly specified in the header.
Bias variations over the night are often modest. We can check whether use of the overscan region is necessary by looking at the statistics of the overscan regions of the bias image, as debiased by the master bias frame alone. If these are near zero with minimal variance, then the master bias frame alone is doing a good job.
I start by reading the keywords from the header that tell me where on in the image the overscan regions are stored, and then extract the corresponding regiens:
>>> dataHdr=pyfits.getheader("mdR-00114174.fits")
>>>
>>> print dataHdr['BIASSEC0']
[21:40,22:2069]
>>> biasOvscn0 = dbImage[21:2069,20:40]
>>>
>>> print dataHdr['BIASSEC1']
[2089:2108,22:2069]
>>> biasOvscn1 = dbImage[21:2069,2088:2108]
Note the difference between the ordering of axis and
pixel indexing and slicing conventions in Python and the FITS header,
whose keywords were written folling IRAF convencions.
(Normally, I would extract the numbers from the strings returned from the header, but the example is easier to follow without the string handling.)
Carefully examine the overscan regions, first by looking at the images:
>>> numdisplay.display(biasOvscn0) Image displayed with Z1: -6.0 Z2: 8.0 >>> numdisplay.display(biasOvscn1) Image displayed with Z1: -2.0 Z2: 10.0Notice that there are signs that something funny is going on; when zoomed out, you can see faint stripeing. Further exploration shows that, although strange, the level of the signal involved is much lower than the noise we expect in our final image, and so we ignore it.
Now we look at some statistics to see whether we need to change the level based on the overscan region, and whether differences in pixel to pixel variations are significant:
>>> biasLevel0 = mean(reshape(biasOvscn0,size(biasOvscn0))) >>> biasStd0 = std(reshape(biasOvscn0,size(biasOvscn0))) >>> expectedStd0 = dataHdr['RDNOISE0']/dataHdr['GAIN0'] >>> print "Amp 0: mean %f, measured std. %f, expected std. %f" % (biasLevel0, biasStd0, expectedStd0) Amp 0: mean 1.954028, measured std. 1.359730, expected std. 1.520000 >>> >>> biasLevel1 = mean(reshape(biasOvscn1,size(biasOvscn1))) >>> biasStd1 = std(reshape(biasOvscn1,size(biasOvscn1))) >>> expectedStd1 = dataHdr['RDNOISE1']/dataHdr['GAIN1'] >>> print "Amp 1: mean %f, measured std. %f, expected std. %f" % (biasLevel1, biasStd1, expectedStd1) Amp 1: mean 4.498950, measured std. 1.476259, expected std. 1.770000 >>>Note that the standard deviation of the bias section is less than than expect from the read noise alone after the master bias is applied; the master bias does a good job taking out the pixel to pixel variations in the bias. So, we only need to apply the offset in bias level:
>>> print dataHdr['DATASEC0'] [41:1064,22:2069] >>> dbImage[:,:1064] = dbImage[:,:1064] - biasLevel0 >>> >>> print dataHdr['DATASEC1'] [1065:2088,22:2069] >>> dbImage[:,1064:] = dbImage[:,1064:] - biasLevel1
We now have our debiased image stored in the dbImage python array.
7.3 Flat fielding
The next step is to "flat field" the image, or rescale each pixel so that the "a" term has a common value for all pixels.
Combine all the flat field exposures and scale the result so that
>>> # Create an empty list in which to accumulate flat field exposures
... flatList = []
>>> # Look up the exposure IDs for the r flat fields from the night log
... flatExpIds = (114044,114049,114054,114060)
>>> for expId in flatExpIds:
... #
... # Read the raw flat field exposure from disk
... #
... flatRaw = pyfits.getdata("mdR-%08d.fits" % expId)
... #
... # The flatfield exposures need to be debiased, just like the data
... # exposure. Here, I assume the overscan regions are the same
... # as in the data exposure. It would be safer to etract the region from
... # each header.
... #
... dbFlat = flatRaw - bias
... biasOvscn0 = dbFlat[21:2069,20:40]
... biasLevel0 = mean(reshape(biasOvscn0,size(biasOvscn0)))
... dbFlat[:,:1064] = dbFlat[:,:1064] - biasLevel0
... biasOvscn1 = dbFlat[21:2069,20:40]
... biasLevel1 = mean(reshape(biasOvscn1,size(biasOvscn1)))
... dbFlat[:,1064:] = dbFlat[:,1064:] - biasLevel1
... #
... # Append our debiased flat field to the list
... #
... flatList.append(dbFlat)
...
>>>
Create a single array that is the median of the flat field exposures:
>>> flat = median(array(flatList))The "array" function turns a list of arrays of the same size and type into a single array, in which the first axis corresponds to the index or the list. The median command takes the median along the first axis.
Now extract just the data section of the flat field. Again, it would be safer to determine the limits using the DATASEC0 and DATASEC1 keywords from the FITS header:
>>> dataFlat = flat[21:2069,40:2088]
Now normalize the flat field to have a mean of 1. The "mean" function only takes the mean along the first axis, so I use the "reshape" function to cast the entire 2-dimensional flat field image into a one dimensional array.
>>> flat = flat / (mean(reshape(dataFlat,size(dataFlat))))
Finally, apply the flat field to data image, and display the result:
>>> ffImage = dbImage/flat >>> numdisplay.display(ffImage,z1=0,z2=1000) Image displayed with Z1: 0 Z2: 1000 >>> numdisplay.display(ffImage,z1=400,z2=500) Image displayed with Z1: 400 Z2: 500
7.4 Write the data image, with necessary changes to the header
We now have the data in our original image bias subtracted and flat fielded. We now need to extract the part of the array which is true sky data:
>>> ffData = ffImage[21:2069,40:2088]
and write it to a FITS file on disk. Most of the metadata in the FITS header in the original exposure applies the the corrected image, so we can copy the header from the original exposure FITS file and then update or remove the invalidated FITS keywords:
>>> ffData = ffImage[21:2069,40:2088]
>>> hdu=pyfits.PrimaryHDU(ffData,dataHdr)
>>> for keyword in ["DATASEC0","DATASEC1","BIASSEC0","BIASSEC1"]:
... del hdu.header[keyword]
...
>>> hdu.header["CRPIX1"] = hdu.header["CRPIX1"] - 40
>>> hdu.header["CRPIX2"] = hdu.header["CRPIX2"] - 21
>>> hdu.header.add_comment("Bias subtracted, flat fielded, and trimmed")
>>> outContent = pyfits.HDUList([hdu])
>>> outContent.writeto("im-114174.fits")
