Tessellating images

PAH feature emission (3290 nm) in the planetary nebula BD+30° 3639

program names, file names
variable names
prompts, commands, program code


Contents:

Overview

Tessellating (or co-adding) the images means lining up all the images, adding them together, and dividing by the number of images. This has the effect of reducing the noise in the background, which enhances the signal-to-noise ratio (SNR). There are at least two factors which contribute to the noise in the images: read noise and shot noise. Shot noise is likely the dominant source of noise, and that will almost certainly be so if the wavelength of the observation is greater than 2 microns. Shot noise ends up in the data because blackbody emission from the atmosphere, the telescope, and the instrument's optics compete with the emission from the astrophysical source being observed. Read noise is a property of the electronics and can sometimes be diminished with larger numbers of Fowler samples in the integrations, but that's something that has to happen at the telescope. So we'll just talk about dealing with shot noise.

The statistics of shot noise is such that the noise is diminished by longer integration times. This is exactly analgous to the result from Gaussian statistics which says that the uncertainty in an experiment decreases with the square root of the number of measurements. In the same way, shot noise decreases with the square root of the total integration time, ie., an image with an integration time of 40 seconds will have approximately twice the SNR of an image with an integration time of 10 seconds. Therefore, when we are observing, we integrate for as long as possible. Ideally, we'd just point the telescope at the object and instruct the camera to integrate for as long as we want. Of course, this doesn't work in practice, principally because the detector eventually saturates (but also because the telescope may not track perfectly, etc.) So instead, we take take a large number of images with shorter integrations. This allows us to accumulate a significant integration time, but it keeps the detector from saturating. The consequence is that we must go through and line up all the images before tessellating them. This is a major pain in the @$$.

Under the best of conditions, the amount of noise in each image is fairly consistent throughout the images we are tessellating, and the calibration of each image is the same. Examples when this isn't true include instances in which it is necessary to tessellate images taken on more than one night, such that the calibration differs from night to night, and perhaps the images taken one night are significantly noisier than images from another night. Another possibility is that the airmass of the object changes enough through the course of the observation as to make a significant difference in the calibration (ie., more than around 10%). Under the best of conditions, the same calibration can be applied to all the images, and the calibration can be applied in a single stroke after tessellating. Sometimes, however, it may be necessary to calibrate each image individually. It may also be helpful to give appropriately-chosen weighting factors to different images based on how noisy the backgrounds are (this is determined, for example, by using the RMS values reported by bkgnorm). As we'll see, both of these special cases can easily be dealt with using the software described below.

As a couple of final points in this overview it's important to make it clear that you should only tessellate images that have the same wavelength. At the end of this part of the data reduction, you'll have one tessellated image for each wavelength, where each tessellated image was created by co-adding several data images at that wavelength. Also, it's generally not necessary to tessellate images of standard stars, as the SNR for those images is usually quite high.


newwin and wblink

In the work we'll be doing next, it's useful to have two display windows on the screen at once, such that the two windows are directly on top of each other, so that we can blink between the two windows. In this way we'll be able to directly compare the alignment of two images.

The startup file (in a graphics session of IDL) will automatically bring up a display window in a particular position on the computer's monitor. This is done by running the newwin procedure, and the code of this program determines where the display window appears on the monitor (feel free to adjust this according to your own preference). Whenever a new display window is created by newwin, the new window becomes the active display window, and any display commands send their output to the active window. You can identify the currently active display window by checking the value of the !d.window environment variable:
IDLprompt> print, !d.window
It is possible, however, to specify the active display window, even if several windows have been created by newwin. The different windows are identified by the number at the top of the window---the original window is likely 32, and successively created windows are probably 33, 34, etc. The following command makes window number 33 the active window:
IDLprompt> wset, 33
(If window number 33 does not currently exist, an error will occur, a message will be printed to the command screen, and nothing else will happen---the active display window will be whatever it was before.)

Another couple of useful procedures are wshow and wdelete. The following command will delete window number 33 (ie., window 33 will vanish from the monitor and no longer be available for displaying images:
IDLprompt> wdelete, 33
(As before, if window number 33 does not currently exist, an error will occur, a message will be printed to the command screen, and nothing else will happen). The following command brings window number 33 to the forefront on the monitor:
IDLprompt> wshow, 33
(And again, if window number 33 does not currently exist, an error will occur, a message will be printed to the command screen, and nothing else will happen). If window number 33 is iconified when you run this command, nothing will happen. The following command will work, even if the window is iconified:
IDLprompt> wshow, 33, iconic=0

If you have two windows currently created, and you want to blink between them, use the wblink procedure. The following command will blink windows 33 and 34:
IDLprompt> wblink, 33, 34
Now we can blink between the two windows by clicking the first two mouse buttons (index finger for number 33 and middle finger for number 34). Clicking the third button (ring finger) ends the wblink program.


Aligning the images

The first thing to do is to determine the registration of all the images to be tessellated (co-added). The registration refers to the process of aligning the images. The easiest way to do this is to use cuts on something bright and compact which appears in each image, preferably a field star (we'll also consider the case where there is no field star a bit later). By running cuts on all the images and getting the coordinates of the field star in each image, we have the two registration parameters (x- and y-coordinates) that we need to align all the images. At this point, you should pick an image which will serve as the anchor image, to which all the other images will be aligned. Ideally, this image will be an image near the middle of the observation (so it serves as sort of an average of the observation conditions, rather than an extreme, like the first or last image in the observation) which has an overall good quality to it, ie., a fairly crisp appearance, with a minimum of blurring due to atmospheric effects, telescope motion, etc. This is largely a judgement call.

Now you should determine the amount by which all the other images should be shifted to line up with the anchor image. This is done, as a first approximation, simply by subtracting the star coordinates of all the images to be shifted from the coordinates of the star in the anchor image:

dXi = Xa - Xi
dYi = Ya - Yi
where Xa is the x-coordinate of the star in the anchor image,Yi is the y-coordinate of the star in the ith image to be shifted, and dXi is the amount by which the ith image should be shifted in the x-direction, etc.

Now that we know approximately how much to shift the images, it's a good idea to make a visual inspection of each image to see how good the alignment is. If the observing conditions were not great, eg., if the field star is blurred, the alignment suggested by cuts may need some adjustment. A straightforward way to accomplish this is to (1) display the anchor image (using v), (2) bring up a new display window in the same place on the monitor (just run newwin with no arguments:
IDLprompt> newwin
and you can do this to create as many windows as you want, within reason), and then (3) display the shifted image in the new window:
IDLprompt> v, shiftxy(impackname, dX, dY, /wt)
shiftxy will shift the image in impackname by the amounts dX and dY, and the shifted image will be displayed (no change is made to the image in impackname). This will work even if dX and dY are not integers (ie., shifting by fraction pixel amounts). Setting the wt keyword applies the weight plane before shifting, which is generally a good idea. Indicentally, shiftxy wraps as it shifts. So if dX=5, for example, then the five columns on the right-hand side of the image will be wrapped to the left-hand side of the displayed image. This wrapping happens only when displaying, not when actually tessellating.

Now we can use wblink to blink between the two images. cuts is generally pretty good at giving alignment parameters if there's a field star, but it may be necessary to adjust dX and dY by eye. Keep adjusting dX and dY and blinking the images until the alignment looks best.

If there is no field star or other fiducial to use for registration purposes, the images will just have to be aligned by eye from the outset. In that case, just pick an image as the anchor, and then use shiftxy and wblink to determine the best shift parameters that you can. It's tedious work, any way you slice it, but it's important to have good shift parameters before trying to tessellate the images.


Co-adding the images

Now it's time to actually co-add the images using the coadd procedure. The first positional parameter can be a 2-D image, a 3-D image (such that the top plane is the image and the bottom plane is the weight plane), or an impack (which itself contains a 2-D or 3-D image). coadd uses weight planes in the process of tessellating, so if coadd can't derive a weight plane from the first input parameter, it uses a plane of ones which is the same size as the image plane.

The next two positional parameters are the shift parameters dX and dY. We'll run coadd once for each image, and we need to set the new keyword for the first, and only the first, image that we co-add. We need to set the neg keyword for any images which are in the background (for observations that were nodded on chip, for example). Finally, we supply an undefined, named variable as an extra positional parameter when we run coadd on the last image to be co-added, and this variable will, after runtime, be a 3-D array consisting of the tessellated image and the tessellated weight plane.

As an example, let's say that we need to co-add the images in four impacks: impack1, impack2, impack3, and impack4, and none of the images were nodded on chip, so all the images are in the foreground. And we'll refer to the shift parameters for the images such that dX1 is the amount by which the image in impack1 should be shifted in the x-direction, etc. We'd use the following four commands to co-add these four images:
IDLprompt> coadd, impack1, dX1, dY1, /new IDLprompt> coadd, impack2, dX2, dY2 IDLprompt> coadd, impack3, dX3, dY3 IDLprompt> coadd, impack4, dX4, dY4, temp
Now temp will have the co-added image and weight plane. You still need to calibrate the image plane, using the calibration factor calculated by jyperadu (we'll pretend that it's 3.5x10-6 Jy/ADU). The easiest way to do this is to put the contents of temp into an impack (let's call it fido), for which we'll need a FITS header (let's say that impack3 was the anchor, so we'll alter its header for our new impack). Then we just use imarith to make the calibration, and wfits to save the calibrated image to the file fido.fts (in the current directory):
IDLprompt> head_edit, filter( impack3 , /h ), newh, add_message = 'coadded' IDLprompt> imcre, fido, filter( temp ) , newh , wt = filter( temp , /wt ) IDLprompt> imarith, fido, 3.5e-6, /multiply, headermessage = 'calibrated' IDLprompt> wfits, 'fido.fts', fido

Now let's consider a more complicated example. Let's say that you have four images to be co-added, but that the observations were made on two nights, with two images for each night. Let's also say that the images were nodded on chip, so that the first two images (from the first night) are in impack1, and the last two images (from the second night) are in impack2. Additionally, lets say that the images from the second night are quite a bit noisier than the images from the first night, which we know from looking at RMS values from the two impacks: RMS1 < RMS2 (recall that the bkgnorm procedure reports RMS values). And because the images were taken on two different nights, there will be two different calibrations: calib1 for the images in impack1 (from the first night), and calib2 for the images in impack2 (from the second night). We'd use the following four commands to co-add these four images:
IDLprompt> coadd, impack1, dX1, dY1, mf=calib1, snr=RMS1, /new IDLprompt> coadd, impack1, dX2, dY2, mf=calib1, snr=RMS1 IDLprompt> coadd, impack2, dX3, dY3, mf=calib2, snr=RMS2 IDLprompt> coadd, impack2, dX4, dY4, temp, mf=calib2, snr=RMS2
Again, temp will have the co-added image and weight plane. This example shows how the snr and mf keywords can be used to apply varying calibration and weighting parameters to the images being tessellated.

This time the image in temp is already calibrated, so we can just make a new impack (after making a header) and write the image to disk:
IDLprompt> head_edit, filter( impack1, /h ), newh, $ IDLprompt> headermessage = 'coadded and calibrated' IDLprompt> imcre, spot, filter( temp ), newh, wt = filter( spot , /wt ) IDLprompt> wfits, 'spot.fts', spot

At this point, we are essentially done with the weight plane. It's a good idea to keep it (as we have done with imcre), but we won't need it again.


Carl Welch