8 - Co-Adding and Dithering

Referenced 06/10/08: http://www.phys.unm.edu/~cpo/html/twhtml/iraf/coadd.html

Using imalign and imcombine for AWO images

Let's say you're doing an imaging experiment, with the goal of achieving some limiting Signal-to-noise or limiting surface brightness. How do you go about obtaining your images? If you calculate that you will need 90 minutes of exposure time at Capilla to achieve these limits, do you

-- take one 90-minute exposure of your object, or
-- take several shorter exposures and combine them?

There are several advantages to doing the latter. In the above example, the strategy might be to take 6-15 minute exposures or maybe 3-30 minute exposures of your object. Between each exposure, you offset the telescope slightly (maybe 10 arcsec) so that on none of the series of images is your object in exactly the same location on the chip. For a sequence of 3 exposures you might offset 10" E of the center image and then 10" W of the center location. For longer sequences you might pattern them in a box or cross-like pattern. In the end, you will simply have to shift all the images until your object lines up at the same pixel location and then combine them using a median-filter coaddition. Why go to this bother?



-- Cosmic rays: by stacking a series of registered images you can reject pixels which have very deviant values from the others in the series (apply a median filter with rejection). Since cosmic rays arrive at random locations on the chip it is very unlikely that two cosmic rays would be found at the exact same location after shifting and registering your images. Thus, bothersome cosmic-ray removal techniques are not necessary.
-- ReadNoise: by median-filtering a set of images you effectively reduce the contribution of the ReadNoise on your images by the square root of N, where N is the number of images in your sequence. This simply results from the fact that for any given pixel, you have sampled it N times in uncorrelated exposures. The error on the mean for a pixel value in this series is then reduced by the square root of N (like measuring a standard deviation).
-- Similarly, you also reduce the uncertainty or noise in the sky-background by the square root of N for the same reasons as above.
-- Errors in flatfielding: depending on your science requirements very accurate flatfielding is often needed. Since you are sampling your object at different locations on the chip, you effectively average out the small-scale flat-fielding errors over the area in which you dithered your exposure.

Do you achieve the same sensitivity with 3 30-min exposures as with 1 90-min exposure? The answer is yes, because the S/N of your image increases with the square root of the exposure time (for sky-noise limited data). Thus, reducing the noise contribution on a series of 3 exposures by root 3 is the same as integrating three times longer. It is important to note that each of the images in your series should be sky-noise dominated rather than read-noise dominated in order to see the gains shown above. The transition from readnoise to sky-noise (Schott-noise) limiting your sensitivity depends on the sky-brightness (phase of moon), the filter (color of the sky), and the amount of readnoise/pixel. You can calculate this using the standard S/N equations.

Here's a typical example of the IRAF procedure necessary to register and combine your images. The tasks are imalign and imcombine in the images.immatch package. Before attempting any of the following steps on your images, be sure that all the steps described earlier in this document have been completed. Image alignment and co-addition should be carried out only on fully reduced images.

Say you have three images of the same object, each shifted a bit from the others...image1, image2, and image3. These image names should be listed in the file "images.list". Select one of the images to be the reference image, the one to which you wish all the other images to align (i.e. centered in the field). Let that image be image1. Also create a coordinate file, using the imexamine program to measure the centroids of 4 or 5 stars in the field on the reference image. Make sure that none of the stars you select are shifted off of the frame because of the dithering. Simply record these centroids (1 per line) in a file called "images.coord". Next, you want to estimate the shifts for the other frames in your sequence. Just choose one of the stars in the images.coord file, and measure it on the other images in the sequence. You estimate the shift by subtracting the location of that star on the reference frame from the other images:

(xref - xobj, yref - yobj). Your lists, contained in 3 different files, should look something like this:

    images.list    images.coord    images.shift
     image1         374  557             0  0
     image2         597  533            64  0
     image3         602  501           -64  0
                    655  453

Note that the first line of the images.shift file indicates no shift on the reference image relative to itself.

 cl > lpar  imalign 
input   =         @images.list    Input images
referenc=         image1          Reference image
coords  =         images.coord    Reference coordinates file
output  =      t//@images.list    Output images
(shifts =         images.shift)   Initial shifts file
(boxsize=                    7)   Size of the footnotesize centering box
(bigbox =                   11)   Size of the big centering box
(negativ=                   no)   Are the features negative ?
(backgro=                INDEF)   Reference background level
(lower  =                INDEF)   Lower threshold for data
(upper  =                INDEF)   Upper threshold for data
(niterat=                    3)   Maximum number of iterations
(toleran=                    0)   Tolerance for convergence
(shiftim=                  yes)   Shift the images ?
(interp=                 poly3)   Interpolant
(boundar=                 wrap)   Boundary type
(constan=                   0.)   Constant for constant boundary extension
(trimima=                   no)   Trim the shifted images ?
(verbose=                  yes)   Print the centers, shifts, and trim section ?
(list   =                     )  
(mode   =                   ql)

The program uses the estimated shifts in the images.shift file as an initial guess in calculating the centroids of all the stars in the coordinate file for each of the images in the sequence. It then calculates the actual shift to use for each frame based on a mean of the measurements. If the errors are more than a few tenths of a pixel, you probably need to pick brighter stars. You can edit the coordinate-file and try again. If you set shiftim=yes the program will output the shifted files prefixed by the letter "t" in the above example. Following this prompt cycle

 im >  imalign
 im >  Input images < (@images.list): 
 im >  Reference image (image1):  im >  Reference coordinates file (images.coord): 
 im >  Output images (t//@images.list):
you'll see something like the following output:

         #Coords    Image      X-center  Err         Y-center Err       Num  
                   image1      376.251 (0.127)       555.754 (0.130)     1  
                   image1      595.936 (0.096)       532.124 (0.101)     2  
                   image1      606.788 (0.097)       504.692 (0.119)     3  
                   image1      660.830 (0.093)       453.433 (0.094)     4  
									    
                   image2      311.840 (0.110)       555.951 (0.125)     1  
                   image2      531.612 (0.126)       531.645 (0.125)     2  
                   image2      542.603 (0.100)       504.540 (0.126)     3  
                   image2      596.583 (0.088)       453.292 (0.112)     4  
									    
                   image3      439.966 (0.142)       555.669 (0.122)     1  
                   image3      660.005 (0.112)       531.517 (0.144)     2  
                   image3      670.753 (0.112)       504.442 (0.175)     3  
                   image3      724.541 (0.111)       453.030 (0.115)     4  


   #Refcoords    Reference    X-center   Err         Y-center  Err       Num 
                   image1      376.251 (0.127)       555.754  (0.130)     1  
                   image1      595.936 (0.096)       532.124  (0.101)     2  
                   image1      606.788 (0.097)       504.692  (0.119)     3  
                   image1      660.830 (0.093)       453.433  (0.094)     4                                                
                                                                             
    #Shifts        Image      X-shift    Err         Y-shift    Err       N      Internal
                   image1       0.00    (0.07)        0.00     (0.08)     4     (0.00,0.00)
                   image2      64.29    (0.07)        0.14     (0.08)     4     (0.05,0.14)
                   image3     -63.86    (0.08)        0.34     (0.09)     4     (0.09,0.11)
                                                                                   

   #Trim-Section = [66:1000,2:1024]

   # Shifting images:   

You then want to combine these registered images into a single frame, say "result", using imcombine:

 
 cl > imcombine  t//@images.list result  comb=median  reject=avsigclip  scale=none  zero=mode  statsec=[x1:x2,y1:y2] 

The "zero=mode" option uses the mode of the region specified in statsec to scale each of the images to the same zero-level or sky background. statsec should be set to a blank region on the images for sky level estimation.

You now have a coadded image of your object.


Using Above Procedure for AWO Observations

Meade LX200GPS + SBIG ST8XE + Equinox 5.3

Using images taken during 05/15/08 of the Ring Nebula

--> cat im.lis
   LightG_NXXX_B2-15C_30S0-1.fits
   LightG_NXXX_B2-15C_30S0-2.fits
   LightG_NXXX_B2-15C_30S0-3.fits
   LightG_NXXX_B2-15C_30S0-4.fits
   LightG_NXXX_B2-15C_30S0-5.fits

--> cat image.coords
445.88 431.98
238.45 472.64
114.57 422.70
63.23 357.71
97.90 137.51
628.96 314.57 --> cat images.shift 0 0 1.01 1.82 2.45 2.44 5.14 2.84 3.24 3.76

The shifts in X and Y pixels were an average for each of the 6 stars used for all 4 images. We're ignoring any rotation of the images and assuming that each star in each image used shifts the same in X and Y. This means your image.list file and your images.shift file will have the same number of lines. Your image.coords file will have as many lines (and corresponding X and Y coordinates) as stars you choose to take part in the alignment program. Use at least 2 stars in different areas of the image.

Parameters for imalign were:

--> lpar imalign
 input 		=@im.lis         Input images
 reference	= LightG_NXXX_B2-15C_30S0-1 Reference image
 coords 	= image.coords    Reference coordinates file
 output 	= t//@im.lis      Output images
 (shifts 	= images.shifts)  Initial shifts file
 (boxsize 	= 7)              Size of the small centering box
 (bigbox 	= 15)             Size of the big centering box
 (negative 	= no)             Are the features negative ?
 (background 	= INDEF)          Reference background level
 (lower 	= INDEF)          Lower threshold for data
 (upper 	= INDEF)          Upper threshold for data
 (niterate 	= 3)              Maximum number of iterations
 (tolerance 	= 0)              Tolerance for convergence
 (maxshift 	= INDEF)          Maximum acceptable pixel shift
 (shiftimages 	= yes)            Shift the images ?
 (interp_type 	= poly3)          Interpolant
 (boundary_type = wrap)           Boundary type
 (constant 	= 0.0)            Constant for constant boundary extension
 (trimimages 	= no)             Trim the shifted images ?
 (verbose 	= yes)            Print the centers, shifts, and trim section ?
 (list = ) 
 (mode = ql) 

Parameters for imcombine were:

--> lpar imcombine
 input		 = t//@im.lis      List of images to combine
 output		 = combined2       List of output images
 (headers	 = )               List of header files (optional)
 (bpmasks	 = )               List of bad pixel masks (optional)
 (rejmasks	 = )               List of rejection masks (optional)
 (nrejmasks	 = )               List of number rejected masks (optional)
 (expmasks	 = )               List of exposure masks (optional)
 (sigmas	 = )               List of sigma images (optional)
 (logfile	 = combined_images) Log file
 (combine	 = median)         Type of combine operation
 (reject	 = avsigclip)      Type of rejection
 (project	 = no)             Project highest dimension of input images?
 (outtype	 = real)           Output image pixel datatype
 (outlimits	 = )               Output limits (x1 x2 y1 y2 ...)
 (offsets	 = none)           Input image offsets
 (masktype	 = none)           Mask type
 (maskvalue	 = 0.0)            Mask value
 (blank		 = 0.0)            Value if there are no pixels
 (scale		 = none)           Image scaling
 (zero		 = mode)           Image zero point offset
 (weight	 = none)           Image weights
 (statsec	 = [350:130,480:220]) Image section for computing statistics
 (expname	 = EXPTIME)        Image header exposure time keyword
 (lthreshold	 = INDEF)          Lower threshold
 (hthreshold	 = INDEF)          Upper threshold
 (nlow		 = 1)              minmax: Number of low pixels to reject
 (nhigh		 = 1)              minmax: Number of high pixels to reject
 (nkeep		 = 1)              Minimum to keep (pos) or maximum to reject (neg)
 (mclip		 = yes)            Use median in sigma clipping algorithms?
 (lsigma	 = 3.0)            Lower sigma clipping factor
 (hsigma	 = 3.0)            Upper sigma clipping factor
 (rdnoise	 = 15)             ccdclip: CCD readout noise (electrons)
 (gain		 = 2.63)           ccdclip: CCD gain (electrons/DN)
 (snoise	 = 0.)             ccdclip: Sensitivity noise (fraction)
 (sigscale	 = 0.1)            Tolerance for sigma clipping scaling corrections
 (pclip		 = -0.5)           pclip: Percentile clipping parameter
 (grow		 = al)