| imcoadd (Oct2004) | gemini.gemtools | imcoadd (Oct2004) |
[0.25*Xmax:0.375*Xmax,0.25*Ymax:0.375*Ymax]
[0.625*Xmax:0.75*Xmax,0.625*Ymax:0.75*Ymax]
where Xmax and Ymax are the sizes of x-axis and the y-axis of the
reference image, respectively.
# trimmed|untrimmed
x1 x2 y1 y2
The first line contains "trimmed" or "untrimmed" depending on
whether the image was trimmed relative to the bad pixels defined in
badpixfile.
If badpixfile="default" the default bad pixel file for the
instrument is used if available. For GMOS-N images, gmos$data/mgmos_bpm11.pl
is used for unbinned data and gmos$data/mgmos_bpm22.pl is used
for data binned 2x2. For NIRI images niri$data/niri_bpm.pl is used.
For NIRI images niri$data/niri_bpm.pl is used.
For Hokupaa+QUIRC images quirc$data/quirc_bpm.pl is used.
If badpixfile="" no bad pixel file is used.
1. General description and guidelines --------------------------------------IMCOADD can be used for
* deriving the geometrical transformation (shift, rotation and scaling) necessary to register several images of the same field * applying the derived transformation to each of the images * deriving an average image of the transformed images, cleaning for cosmic-ray-eventsIn order for IMCOADD to work best the images need to contain 5 or more point-like sources. It will in general not work on star forming regions where most of the field is filled with diffuse emission from gas and dust.
IMCOADD is modular in design and each step has a flag in the parameter file related to it. The steps and flags are
1. Find the objects in the reference image (the first in the list of
images) using DAOFIND: fl_find+
2. Derive the transformations using GEOMAP: fl_map+
3. Transform the images using GEOTRAN: fl_trn+
4. Derive the median image using IMCOMBINE: fl_med+
5. Derive the average image cleaned for cosmic-ray-events, using
IMCOMBINE: fl_add+
6. Derive the average uncleaned image using IMCOMBINE: fl_avg+
The steps should be run in order of the flags. If IMCOADD is used with alignmethod=user the most efficient way to run the full process is to call the script twice as follows.
1. First call with the flags set to
fl_find+ fl_map+ fl_trn- fl_med- fl_add- fl_avg- [fl_inter+]
2. Second call with the flags set to
fl_find- fl_map- fl_trn+ fl_med+ fl_add+ fl_avg+
If the input images contain few objects and DAOFIND cannot reliable be made to find them, the positions can be put in a file in the format
x y NoError object-number
(x,y) are the pixel coordinates, "NoError" makes IMCOADD act as if the coordinate file was from DAOFIND, and "object-number" is an object number running from 1 to the total number of objects in the image. The file should be named [reference]_pos, where "reference" is the name of the first image in the list of images to process.
Then run IMCOADD with the flags set to
fl_find- fl_map+ fl_trn+ fl_med+ fl_add+ fl_avg+ [fl_inter+]
IMCOADD can be started from any intermediate result by setting the preceding flags to "no" and the flags for the following actions to "yes". For example, if to rederive the cleaned average image rerun the task with the flags set to
fl_find- fl_map- fl_trn- fl_med- fl_add+ fl_avg+ fl_overwrite+
fl_overwrite=yes will delete the existing output files that are about the be rederived.
2. Input files ---------------The input images should be given as a comma separated list, e.g.
images=image1,image2,image3,image4
The input images can also be given as a file that lists one image per
line, e.g.
images=@imagelist
where "imagelist" contains the image names.
If the input images are MEF files and contain data quality planes, these can be used as the bad pixel masks by setting immasks="DQ". This is the default value of the parameter. If the input images are not MEF files, immasks="DQ" is ignored by the task. Individual bad pixel files can be given for the images to facilitate masking of pixels that are only bad in one of the images. If individual bad pixel files are used, the number of entries in the parameter immasks need to be the same as the number of input images, e.g.
immasks="image1.msk,image2.msk,image3.msk,image4.msk"
If not all the images require an individual bad pixel file, any of the
entries can be substituted with the word "none",
e.g.
immasks="image1.msk,image2.msk,none,image4.msk"
The individual bad pixel files can also be given as a file that lists
one bad pixel file per line, e.g.
immasks=@masklist
where "masklist" contains the bad pixel file names.
The individual bad pixel files may either be IRAF mask files, or ascii files. IRAF mask files are assumed to have "0" in the good pixels and values larger than zero in the bad pixels. If ascii files are used for the individual bad pixel files, then the format must be as show in this example which masks the image sections [1:50,1:50] and the pixel (250,235)
# trimmed
1 50 1 50
250 250 235 235
If the first line in the file contains "# trimmed", IMCOADD will use the file as it is. If the first line in the file contains "# untrimmed", imcoadd will use the image header keyword TRIMSEC to correct the mask file according to how the image was trimmed. In this case, make sure that your images still contain the header keyword TRIMSEC after the trimming.
A bad pixel file for the detector may be given either as an ASCII file in the same format as used for the individual bad pixel files or as an IRAF mask file.
If you run with fl_find- fl_map+, IMCOADD will look for the coordinate file [reference]_pos, where "reference" is the first image in the list of input images.
3. Output files ----------------Finding the objects (fl_find+) will result in a coordinate file named [reference]_pos.
Deriving the transformations (fl_map+) will result in a database file database. The file is ASCII, and if necessary any failed fits may be deleted from it (see GEOMAP for details). Other output files from this step are
[image]_pos Coordinate files for the images, approximate positions
[reference]_cen Coordinate file for the reference image, positions
derived using APPHOT.CENTER with [reference]_pos as
the input file.
[image]_cen Coordinate files for the other images, positions derived
using APPHOT.CENTER with [image]_pos as the input file.
[image]_trn Position input file for GEOMAP. The entry in database
will also be called [image]_trn. The file is made from
[reference]_cen and [image]_cen
Applying the transformations (fl_trn+) results in transformed images named [image]_trn. These images are sky subtracted. In addition, the sky subtracted reference image [reference]_trn is derived. The subtracted sky level is stored in the header keyword MED_SKY in [image]_trn and [reference]_trn. If the input images are MEF, the transformed images [image]_trn and [reference]_trn will be simple FITS.
Deriving the median image (fl_med+) results in the following output files.
[image]_trn_mag Photometry information used for scaling of the
intensity of the images.
[reference]_mag.tab All the photometry information used for scaling
the intensity of the images.
[reference]_med Median image. If the input images are MEF this image
will be simple FITS.
Deriving the average image cleaned for cosmic-ray-events (fl_add+) results in the following output files.
[reference]badpix.pl Bad pixel mask for reference image
[image]badpix.pl Bad pixel mask for other images
[reference]_add Cosmic-ray cleaned average image. If the input
images are MEF this image will also be MEF.
If outimage!="" [reference]_add will be
renamed to outimage.
Deriving the uncleaned average image (fl_avg+) results in the output image [reference]_avg. If the input images are MEF this image this image will be simple FITS.
The median image should not be used for photometry. It should be considered an intermediate step on the way to the cosmic-ray cleaned average image. The uncleaned average image should only be used for checking the output from the cosmic-ray cleaning.
Note, that if the input images are MEF files, the only output image that is also MEF is [reference]_add. This is also the only image that contains complete header information.
4. Detailed description of the method --------------------------------------4.1. Find the objects: fl_find+
IMCOADD uses DAOFIND for finding the objects in the reference image. IMCOADD will determine the sky level using "midpt" from IMSTATISTICS. The noise in the background is derived from the read-noise, gain and the sky level. datamin and datamax are used as input for DAOFIND. datamin should be smaller than zero, to ensure correct behavior of later steps.
4.2. Derive the transformations: fl_map+
IMCOADD uses GEOMAP for deriving the transformations that will register the images. These are the steps if alignmethod=user:
1. IMCOADD centers the objects in the reference image.
2. IMCOADD displays the reference image.
3. The user points to one (two if rotate=yes or scale=yes) object(s)
to be used as common object(s) for all the images.
4. IMCOADD displays the next image, called "current image" in the following.
5. The user points to the same one/two object(s) in common.
6. IMCOADD makes a first simple fit of the transformation (shift,
shift+rotation, or shift+rotation+scaling).
7. IMCOADD matches the objects in [reference]_pos to the positions in
the current image.
8. IMCOADD derives accurate coordinates for the objects in the current
image using APPHOT.CENTER.
9. IMCOADD uses GEOMAP iteratively to derive the transformation, using all
the objects in common that were successfully centered in both images
The parameters geofitgeom, order, sigfit, niter,
and coolimit are used in this process.
10. If fl_inter=yes, the user is after the last iteration left with
geomap in interactive mode and can inspect the fit, change the fitting
parameters and refit if necessary. A short help menu is printed.
When happy, the users exits geomap with "q".
11. The steps 4-10 are repeated for the rest of the images.
If alignmethod=twodx, the steps 2-6 are replaced by cross-correlation of two image sections specified in asection. The cross-correlation is done with IMMATCH.XREGISTER. The output from the cross-correlation is used for the first simple fit in step 6. It is important to ensure that xwindow is at least twice the maximum pixel offsets between the successive input images. Otherwise XREGISTER will fail to find the cross-correlation peak. Further xwindow cannot be larger than the minimum size of the image sections defined in asection. xwindow will be automatically reset to a smaller value if the input value is too large.
If alignmethod=header, the steps 2-6 are replaced by reading the relevant header keywords and deriving the first rough estimates of the shifts in the X- and Y-directions from the values of these keywords.
If alignmethod=wcs, the steps 2-6 are replaced by deriving the transformations based on the WCS in the headers.
If niter=0, the transformation will not be improved past the first rough estimate. This may be useful if the individual frames do not contain any visible objects, but the header information is accurate enough to allow co-addition of the images.
4.3 Transform the images: fl_trn+
IMCOADD uses GEOTRAN to transform the images according to the derived transformations. Large values of geonxblock and geonyblock causes the task to use more memory and run faster.
If geointer="nearest" no interpolation of pixel values is done. The noise characteristics of the transformed images is the same as for the input images and no correlated noise between the pixel values result from the transformation.
IMCOADD derives the sky level of each of the images using IMSTATISTICS iteratively. The resulting sky level is the "midpt" from IMSTATISTICS for values within 5 times the expected noise in the background. The expected noise is determined from the sky level, the read-noise and the gain. The sky level is subtracted from each of the images. The images are still called [reference]_trn and [image]_trn. The subtracted sky level is recorded in the headers in the keyword MED_SKY.
4.4 Derive the median image: fl_med+
IMCOADD derives the relative signal in each of the images from aperture photometry of the objects that were successfully centered in all the images and that do not contain pixels outside the range defined by datamin and datamax. datamin should be smaller than zero. datamax should be the saturation level of the detector. The aperture size should be large enough to eliminate differences in the photometry due to small differences in the seeing, but small enough to ensure that sky noise does not dominate (as it may in the near-IR).
It is recommended to scale the images in intensity, fl_scale=yes. If fl_scale=no is used, even small differences in the signal between images will show up in the combined image as "missing" or "extra" signal in objects where some pixels in some of the input images are flagged as bad (either in the detector bad pixel file or as a cosmic-ray hit). To illustrate this, imagine a set of images where one image has lower signal that the rest. An extended object with a few bad pixels in the low-signal image will then on the co-added image have slightly higher signal in that area (where the result is the mean of all the high signals) than in the surrounding area (where the result is the mean of all the high signals and one lower signal).
IMCOADD derives the median image, with the images scaled according to their relative signal if fl_scale=yes. The bad pixels as defined in [image]badpix.pl are excluded from the median image.
IMCOADD will automatically turn off the scaling (setting fl_scale=no) if no objects in common are found in the images.
4.5 Derive the cleaned average image: fl_add+
For each input image are calculated
* A difference image: "difference image" = "transformed input image" - "median image" * A sigma image representing the local noise in the difference images
Scalings and sky levels are taken into account. The ron and the gain are needed for these calculations. If scalenoise > 0, then a noise term of "input image"*scalenoise is added in quadrature to the sigma image. The purpose of scalenoise is to make it possible to combine images with significant differences in the PSFs. scalenoise=0 works for images with very little difference in the FWHM of the PSFs, while scalenoise=0.3 can handle images with up to 15-20% difference in the FWHM of the PSFs. When combining images with 15-20% difference in the FWHM of the PSFs using scalenoise=0 the result may be that rings around the objects in the image with the worst seeing get flagged as bad pixels.
Based on difference and the sigma images the deviating pixels are identified. Deviating pixels satisfy one of the following criteria where N(image) is the number of ADU in a given pixel in an image:
N(median) <= limit && N(difference) >= lowsigma * N(sigma)
N(median) <= limit && N(difference) >= lowlimit
For median signal below limit pixels are flagged as deviating
if they deviate more than lowsigma times the local noise, or
with more than an absolute value of lowlimit.
For median signal above limit no flagging is done.
limit should be set to approximately 15 times the noise in the
background. If the images have very different sky levels, it is
recommended that to use an absolute limit set to 15 times the
noise in the sky of the image with the lowest sky level.
In that case set key_limit=no.
If growthrad=1 then pixels neighboring the deviating pixels are also identified as deviating.
The mask image [image]badpix.pl for each input image is updated to include the deviating pixels as bad pixels.
IMCOADD derives the average image, with the images scaled according to their relative signal if fl_scale=yes. IMCOADD will automatically turn off the scaling (setting fl_scale=no) if no objects in common are found in the images.
The bad pixels as defined in [image]badpix.pl are excluded from the average image. It is recommended to inspect the output [image]badpix.pl for possible problems related to the cleaning of the images. If there are rings of excluded pixels around all the stars in the images, then the difference in FWHM combined with the choice of limit made IMCOADD exclude pixels from the images with the worst seeing. The solution to this is to use scalenoise=0.3 and/or a lower value of limit.
4.6 Derive the uncleaned average image: fl_avg+
IMCOADD derives the average image, with the images scaled according to their relative signal if fl_scale=yes. IMCOADD will automatically turn off the scaling (setting fl_scale=no) if no objects in common are found in the images.
No pixels are excluded, and no cleaning takes place.
5. Photometry from the combined cleaned image ----------------------------------------------The median image should not be used for photometry. It should be considered an intermediate step on the way to the cosmic-ray cleaned average image. The uncleaned average image should only be used for checking the output from the cosmic-ray cleaning.
Science grade photometry should be derived from the cleaned average image [reference]_add. [reference]_add is on the same magnitude scale as the reference image, e.g. same exposure time and same airmass. If the reference image was obtained under photometric conditions, then photometry derived from [reference]_add can be calibrated in the same way as would be done for photometry from the reference image.
To improve the absolute accuracy of the photometry derived from [reference]_add, the user may also derive relative offsets between [reference]_add and the input images. This is easiest done by deriving the photometry from [image]_trn, since these are all registered and therefore the same object coordinate file may be used for all the images. A planned upgrade to IMCOADD includes deriving all the magnitude offsets between [reference]_add and the input images.
6. Planned upgrades to IMCOADD ------------------------------- The list is in the order of priority.1. Optional use of the variance plane in MEF files. 2. Output files that are larger than the input files, instead of being cropped to the size of the reference image. 3. Other algorithms for the rough alignment of the images, e.g. use of a user-supplied list of offsets
7. Reference for the procedure use in IMCOADD -----------------------------------------------The procedure implemented in IMCOADD is described in
Jorgensen I., 2003, AJ, submitted
1. Combine the images in "imagelist" using default values of all parameters
cl> imcoadd @imagelist
Hardware: Ultra 30, CPU 300 Mhz, 1Gb RAM at 100 Mhz,
SCSI3 hard drive
Combination of five 1kx1k NIRI images, immasks="", alignmethod=wcs and verbose=yes. Time for completion of all steps: 2 min 43 sec.
Combination of eight GMOS images binned 2x2, immasks="", alignmethod=wcs and verbose=yes. Time for completion of all steps: 23 min.
IMCOADD works on IRAF images and simple FITS. It also works on MEF files that have one science extension. Processing of MEF files requires imtype="fits".
If IMCOADD crashes a number of temporary files and images may be left in the current directory. These have names starting with "tmp" and should all be deleted before attempting another run.
IMCOADD crashes when any of the objects appear to have negative flux. This may happen if the objects are very faint and the sky level gets determined slightly too high. This bug can be avoided by setting the threshold high enough. For manually provided positions for very faint objects, this bug may be difficult to avoid and difficult to fix, since the reason is an untrapped floating point exception in IRAF.
Image names containing a '.' other than for the file name extension should be avoided as in some cases IMCOADD does not parse such image name correctly.