combine: Combine CCD images

Package: ccdred

Usage

combine input output

Parameters

input: List of CCD images to combine. Images of a particular CCD image type may be selected with the parameter ccdtype with the remaining images ignored.

output: Output combined image or list of images. If the project parameter is no (the typical case for CCD acquisition) then there will be one output image or, if the subsets parameter is selected, one output image per subset. If the images consist of stacks then the project option allows combining each input stack into separate output images as given by the image list.

plfile = "" (optional): Output pixel list file or list of files. If no name is given or the list ends prematurely then no file is produced. The pixel list file is a map of the number of pixels rejected or, equivalently, the total number of input images minus the number of pixels actually used. The file name is also added to the output image header under the keyword BPM.

sigma = "" (optional): Output sigma image or list of images. If no name is given or the list ends prematurely then no image is produced. The sigma is standard deviation, corrected for a finite population, of the input pixel values (excluding rejected pixels) about the output combined pixel values.

ccdtype = "": CCD image type to combine. If specified only input images of the specified type are combined. See ccdtypes for the possible image types.

amps = yes: Combine images by amplifier? If yes then the input images are grouped by the amplifier parameter and each group combined into a separate output image. The amplifier identifier is appended to the output image name(s). See subsets for more on the amplifier parameter.

subsets = no: Combine images by subset parameter? If yes then the input images are grouped by subset parameter and each group combined into a separate output image. The subset identifier is appended to the output image name(s). See subsets for more on the subset parameter.

delete = no: Delete input images after combining? Only those images combined are deleted.

clobber = no: Clobber existing output images? THIS OPTION IS NO LONGER SUPPORTED BUT THE PARAMETER REMAINS FOR NOW FOR BACKWARD COMPATIBILITY. IF SET TO yes AN ERROR ABORT WILL OCCUR.

combine = "average" (average|median): Type of combining operation performed on the final set of pixels (after offsetting, masking, thresholding, and rejection). The choices are "average" or "median". The median uses the average of the two central values when the number of pixels is even.

Type of rejection operation performed on the pixels remaining after offsetting, masking and thresholding. The algorithms are discussed in the DESCRIPTION section. The rejection choices are:

     none - No rejection
   minmax - Reject the nlow and nhigh pixels
  ccdclip - Reject pixels using CCD noise parameters
 crreject - Reject only positive pixels using CCD noise parameters
  sigclip - Reject pixels using a sigma clipping algorithm
avsigclip - Reject pixels using an averaged sigma clipping algorithm
    pclip - Reject pixels using sigma based on percentiles

project = no: Project (combine) across the highest dimension of the input images? If no then all the input images are combined to a single output image. If yes then the highest dimension elements of each input image are combined to an output image and optional pixel list and sigma images. Each element of the highest dimension may have a separate offset but there can only be one mask image.

outtype = "real" (short|ushort|integer|long|real|double): Output image pixel datatype. The pixel datatypes are "double", "real", "long", "integer", unsigned short ("ushort") and "short" with highest precedence first. If none is specified then the highest precedence datatype of the input images is used. A mixture of short and unsigned short images has a highest precedence of integer. The datatypes may be abbreviated to a single character.

offsets = "none" (none|wcs|grid|<filename>)

Integer offsets to add to each image axes. The options are:

"none": No offsets are applied.

"wcs": The world coordinate system (wcs) in the image is used to derive the offsets. The nearest integer offset that matches the world coordinate at the center of the first input image is used.

"grid"

A uniform grid of offsets is specified by a string of the form

grid [n1] [s1] [n2] [s2] ...

where ni is the number of images in dimension i and si is the step in dimension i. For example "grid 5 100 5 100" specifies a 5x5 grid with origins offset by 100 pixels.

<filename>: The offsets are given in the specified file. The file consists of one line per image with the offsets in each dimension forming the columns.

masktype = "none" (none|goodvalue|badvalue|goodbits|badbits): Type of pixel masking to use. If "none" then no pixel masking is done even if an image has an associated pixel mask. The other choices are to select the value in the pixel mask to be treated as good (goodvalue) or bad (badvalue) or the bits (specified as a value) to be treated as good (goodbits) or bad (badbits). The pixel mask file name comes from the image header keyword BPM. Note that when combining images by projection of the highest dimension only one pixel mask is applied to all the images. Also if the number of input images becomes too large (currently about 115 .imh or 57 .hhh images) then the images are temporarily stacked and combined by projection which also means the bad pixel mask from the first image will be used for all images.

maskvalue = 0: Mask value used with the masktype parameter. If the mask type selects good or bad bits the value may be specified using IRAF notation for decimal, octal, or hexadecimal; i.e 12, 14b, 0cx to select bits 3 and 4.

blank = 0.: Output value to be used when there are no pixels.

scale = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>): Multiplicative image scaling to be applied. The choices are none, scale by the mode, median, or mean of the specified statistics section, scale by the exposure time in the image header, scale by the values in a specified file, or scale by a specified image header keyword. When specified in a file the scales must be one per line in the order of the input images.

zero = "none" (none|mode|median|mean|@<file>|!<keyword>): Additive zero level image shifts to be applied. The choices are none or shift by the mode, median, or mean of the specified statistics section, shift by values given in a file, or shift by values given by an image header keyword. When specified in a file the zero values must be one per line in the order of the input images. File or keyword zero offset values do not allow a correction to the weights.

weight = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>): Weights to be applied during the final averaging. The choices are none, the mode, median, or mean of the specified statistics section, the exposure time, values given in a file, or values given by an image header keyword. When specified in a file the weights must be one per line in the order of the input images and the only adjustment made by the task is for the number of images previously combined. In this case the weights should be those appropriate for the scaled images which would normally be the inverse of the variance in the scaled image.

statsec = "": Section of images to use in computing image statistics for scaling and weighting. If no section is given then the entire region of the input is sampled (for efficiency the images are sampled if they are big enough). When the images are offset relative to each other one can precede the image section with one of the modifiers "input", "output", "overlap". The first interprets the section relative to the input image (which is equivalent to not specifying a modifier), the second interprets the section relative to the output image, and the last selects the common overlap and any following section is ignored.

Algorithm Parameters

lthreshold = INDEF, hthreshold = INDEF: Low and high thresholds to be applied to the input pixels. This is done before any scaling, rejection, and combining. If INDEF the thresholds are not used.

nlow = 1, nhigh = 1 (minmax): The number of low and high pixels to be rejected by the "minmax" algorithm. These numbers are converted to fractions of the total number of input images so that if no rejections have taken place the specified number of pixels are rejected while if pixels have been rejected by masking, thresholding, or nonoverlap, then the fraction of the remaining pixels, truncated to an integer, is used.

nkeep = 1: The minimum number of pixels to retain or the maximum number to reject when using the clipping algorithms (ccdclip, crreject, sigclip, avsigclip, or pclip). When given as a positive value this is the minimum number to keep. When given as a negative value the absolute value is the maximum number to reject. If there are fewer pixels at some point due to offsetting, thresholding, or masking then if the number to keep (positive nkeep) is greater than the number of pixels no pixels will be rejected and if the number to reject is given (negative nkeep) then up to that number may be rejected.

mclip = yes (ccdclip, crreject, sigclip, avsigcliip): Use the median as the estimate for the true intensity rather than the average with high and low values excluded in the "ccdclip", "crreject", "sigclip", and "avsigclip" algorithms? The median is a better estimator in the presence of data which one wants to reject than the average. However, computing the median is slower than the average.

lsigma = 3., hsigma = 3. (ccdclip, crreject, sigclip, avsigclip, pclip): Low and high sigma clipping factors for the "ccdclip", "crreject", "sigclip", "avsigclip", and "pclip" algorithms. They multiply a "sigma" factor produced by the algorithm to select a point below and above the average or median value for rejecting pixels. The lower sigma is ignored for the "crreject" algorithm.

rdnoise = "0.", gain = "1.", snoise = "0." (ccdclip, crreject)

CCD readout noise in electrons, gain in electrons/DN, and sensitivity noise as a fraction. These parameters are used with the "ccdclip" and "crreject" algorithms. The values may be either numeric or an image header keyword which contains the value. The noise model for a pixel is:

variance in DN = (rdnoise/gain)^2 + DN/gain + (snoise*DN)^2
variance in e- = (rdnoise)^2 + (gain*DN) + (snoise*(gain*DN))^2
               = rdnoise^2 + Ne + (snoise * Ne)^2

where DN is the data number and Ne is the number of electrons. Sensitivity noise typically comes from noise introduced during flat fielding.

sigscale = 0.1 (ccdclip, crreject, sigclip, avsigclip): This parameter determines when poisson corrections are made to the computation of a sigma for images with different scale factors. If all relative scales are within this value of unity and all relative zero level offsets are within this fraction of the mean then no correction is made. The idea is that if the images are all similarly though not identically scaled, the extra computations involved in making poisson corrections for variations in the sigmas can be skipped. A value of zero will apply the corrections except in the case of equal images and a large value can be used if the sigmas of pixels in the images are independent of scale and zero level.

pclip = -0.5 (pclip): Percentile clipping algorithm parameter. If greater than one in absolute value then it specifies a number of pixels above or below the median to use for computing the clipping sigma. If less than one in absolute value then it specifies the fraction of the pixels above or below the median to use. A positive value selects a point above the median and a negative value selects a point below the median. The default of -0.5 selects approximately the quartile point. See the DESCRIPTION section for further details.

grow = 0: Number of pixels to either side of a rejected pixel along image lines to also be rejected. This applies only to pixels rejected by one of the rejection algorithms and not the masked or threshold rejected pixels.

PACKAGE PARAMETERS

The package parameters are used to specify verbose and log output and the instrument and header definitions.

Description

A set of CCD images are combined by weighted averaging or medianing. Pixels may be rejected from the combining by using pixel masks, threshold levels, and rejection algorithms. The images may be scaled multiplicatively or additively based on image statistics, image header keywords, or text files before rejection. The images may be combined with integer pixel coordinate offsets to produce an image bigger than any of the input images. This task is a variant of the images.imcombine task specialized for CCD images.

The input images to be combined are specified by a list. A subset or subsets of the input list may be selected using the parameters ccdtype and subsets. The ccdtype parameter selects only images of a specified standard CCD image type. The subsets parameter breaks up the input list into sublists of common subset parameter (filter, grating, etc.). For more information see ccdtypes and subsets. This selection process is useful with wildcard templates to combine, for example, the flat field images for each filter in one step (see flatcombine). When subsets of the input list are used the output image and optional pixel file and sigma image are given by root names with an amplifier and subset identifier appended by the task.

If the project parameter is yes then the highest dimension elements of each input image are combined to make an output image of one lower dimension. There is no limit to the number of elements combined in this case. This case is If the project is no then the entire input list is combined to form a single output image per subset. In this case the images must all have the same dimensionality but they may have different sizes. There is a software limit of approximately 100 images in this case.

The output image header is a copy of the first image in the combined set. In addition, the number of images combined is recorded under the keyword NCOMBINE, the exposure time is updated as the weighted average of the input exposure times, and any pixel list file created is recorded under the keyword BPM. The output pixel type is set by the parameter outtype. If left blank then the input datatype of highest precision is used. A mixture of short and unsigned short images has a highest precision of integer.

In addition to one or more output combined images there may also be a pixel list image containing the number of pixels rejected at each point in the output image, an image containing the sigmas of the pixels combined about the final output combined pixels, and a log file. The pixel list image is in the compact pixel list format which can be used as an image in other programs. The sigma computation is the standard deviation corrected for a finite population (the n/(n-1) factor) including weights if a weighted average is used.

Other input/output parameters are delete and clobber. The delete parameter may be set to "yes" to delete the input images used in producing an output image after it has been created. This is useful for minimizing disk space, particularly with large sets of calibration images needed to achieve high statistical accuracy in the final calibration image. The clobber parameter allows the output image names to be existing images which are overwritten (at the end of the operation).

An outline of the steps taken by the program is given below and the following sections elaborate on the steps.

o   Set the input image offsets and the final output image size.
o   Set the input image scales and weights
o   Write the log file output

For each output image line:

o   Get input image lines that overlap the output image line
o   Reject masked pixels
o   Reject pixels outside the threshold limits
o   Reject pixels using the specified algorithm
o   Reject neighboring pixels along each line
o   Combine remaining pixels using the weighted average or median
o   Compute sigmas of remaining pixels about the combined values
o   Write the output image line, rejected pixel list, and sigmas

OFFSETS

The images to be combined need not be of the same size or overlap. They do have to have the same dimensionality which will also be the dimensionality of the output image. Any dimensional images supported by IRAF may be used. Note that if the project flag is yes then the input images are the elements of the highest dimension; for example the planes of a three dimensional image.

The overlap of the images is determined by a set of integer pixel offsets with an offset for each dimension of each input image. For example offsets of 0, 10, and 20 in the first dimension of three images will result in combining the three images with only the first image in the first 10 colums, the first two images in the next 10 columns and all three images starting in the 31st column. At the 31st output column the 31st column of the first image will be combined with the 21st column of the second image and the 1st column of the third image.

The output image size is set by the maximum extent in each dimension of any input image after applying the offsets. In the above example if all the images have 100 columns then the output image will have 130 columns corresponding to the 30 column offset in the third image.

The input image offsets are set using the offset parameter. There are four ways to specify the offsets. If the word "none" or the empty string "" are used then all offsets will be zero and all pixels with the same coordinates will be combined. The output image size will be equal to the biggest dimensions of the input images.

If "wcs" offsets are specified then the world coordinate systems (wcs) in the image headers are used to derived the offsets. The world coordinate at the center of the first input image is evaluated. Then integer pixel offsets are determined for each image to bring the same world coordinate to the same point. Note the following caveats. The world coordinate systems must be of the same type, orientation, and scale and only the nearest integer shift is used.

If the input images have offsets in a regular grid or one wants to make an output image in which the input images are "mosaiced" together in a grid then the special offset string beginning with the word "grid" is used. The format is

grid [n1] [s1] [n2] [s2] ...

where ni is the number of images in dimension i and si is the step in dimension i. For example "grid 5 100 5 100" specifies a 5x5 grid with origins offset by 100 pixels. Note that one must insure that the input images are specified in the correct order. This may best be accomplished using a "@" list. One useful application of the grid is to make a nonoverlapping mosaic of a number of images for display purposes. Suppose there are 16 images which are 100x100. The offset string "grid 4 101 4 101" will produce a mosaic with a one pixel border having the value set by blank parameter between the images.

The offsets may be defined in a file by specifying the file name in the offset parameter. (Note that the special file name STDIN may be used to type in the values terminated by the end-of-file character). The file consists of a line for each input image. The lines must be in the same order as the input images and so an "@" list may be useful. The lines consist of whitespace separated offsets one for each dimension of the images. In the first example cited above the offset file might contain:

0 0
10 0
20 0

where we assume the second dimension has zero offsets.

The offsets need not have zero for one of the images. The offsets may include negative values or refer to some arbitrary common point. When the offsets are read by the program it will find the minimum value in each dimension and subtract it from all the other offsets in that dimension. The above example could also be specified as:

225 15
235 15
245 15

There may be cases where one doesn't want the minimum offsets reset to zero. If all the offsets are positive and the comment "# Absolute" appears in the offset file then the images will be combined with blank values between the first output pixel and the first overlapping input pixel. Continuing with the above example, the file

# Absolute
10 10
20 10
30 10

will have the first pixel of the first image in the 11th pixel of the output image. Note that there is no way to "pad" the other side of the output image.

SCALES AND WEIGHTS

In order to combine images with rejection of pixels based on deviations from some average or median they must be scaled to a common level. There are two types of scaling available, a multiplicative intensity scale and an additive zero point shift. The intensity scaling is defined by the scale parameter and the zero point shift by the zero parameter. These parameters may take the values "none" for no scaling, "mode", "median", or "mean" to scale by statistics of the image pixels, "exposure" (for intensity scaling only) to scale by the exposure time keyword in the image header, any other image header keyword specified by the keyword name prefixed by the character '!', and the name of a file containing the scale factors for the input image prefixed by the character '@'.

Examples of the possible parameter values are shown below where "myval" is the name of an image header keyword and "scales.dat" is a text file containing a list of scale factors.

scale = none            No scaling
zero = mean             Intensity offset by the mean
scale = exposure        Scale by the exposure time
zero = !myval           Intensity offset by an image keyword
scale = @scales.dat     Scales specified in a file

The image statistics factors are computed by sampling a uniform grid of points with the smallest grid step that yields less than 10000 pixels; sampling is used to reduce the time need to compute the statistics. If one wants to restrict the sampling to a region of the image the statsec parameter is used. This parameter has the following syntax:

[input|output|overlap] [image section]

The initial modifier defaults to "input" if absent. The modifiers are useful if the input images have offsets. In that case "input" specifies that the image section refers to each input image, "output" specifies that the image section refers to the output image coordinates, and "overlap" specifies the mutually overlapping region of the input images. In the latter case an image section is ignored.

The statistics are as indicated by their names. In particular, the mode is a true mode using a bin size which is a fraction of the range of the pixels and is not based on a relationship between the mode, median, and mean. Also masked pixels are excluded from the computations as well as during the rejection and combining operations.

The "exposure" option in the intensity scaling uses the exposure time from the image header. If one wants to use a nonexposure time image header keyword the !<keyword> syntax is available.

If both an intensity scaling and zero point shift are selected the multiplicative scaling is done first. Use of both makes sense if the intensity scaling is the exposure time to correct for different exposure times and then the zero point shift allows for sky brightness changes.

The image statistics and scale factors are recorded in the log file unless they are all equal, which is equivalent to no scaling. The intensity scale factors are normalized to a unit mean and the zero point shifts are adjust to a zero mean. When the factors are specified in an @file or by a keyword they are not normalized.

Scaling affects not only the mean values between images but also the relative pixel uncertainties. For example scaling an image by a factor of 0.5 will reduce the effective noise sigma of the image at each pixel by the square root of 0.5. Changes in the zero point also changes the noise sigma if the image noise characteristics are Poissonian. In the various rejection algorithms based on identifying a noise sigma and clipping large deviations relative to the scaled median or mean, one may need to account for the scaling induced changes in the image noise characteristics.

In those algorithms it is possible to eliminate the "sigma correction" while still using scaling. The reasons this might be desirable are 1) if the scalings are similar the corrections in computing the mean or median are important but the sigma corrections may not be important and 2) the image statistics may not be Poissonian, either inherently or because the images have been processed in some way that changes the statistics. In the first case because computing square roots and making corrections to every pixel during the iterative rejection operation may be a significant computational speed limit the parameter sigscale selects how dissimilar the scalings must be to require the sigma corrections. This parameter is a fractional deviation which, since the scale factors are normalized to unity, is the actual minimum deviation in the scale factors. For the zero point shifts the shifts are normalized by the mean shift before adjusting the shifts to a zero mean. To always use sigma scaling corrections the parameter is set to zero and to eliminate the correction in all cases it is set to a very large number.

If the final combining operation is "average" then the images may be weighted during the averaging. The weights are specified in the same way as the scale factors. In addition the NCOMBINE keyword, if present, will be used in the weights. The weights, scaled to a unit sum, are printed in the log output.

The weights are only used for the final weighted average and sigma image output. They are not used to form averages in the various rejection algorithms. For weights in the case of no scaling or only multiplicative scaling the weights are used as given or determined so that images with lower signal levels will have lower weights. However, for cases in which zero level scaling is used and the zero levels are determined from image statistics (not from an input file or keyword) the weights are computed from the initial weights (the exposure time, image statistics, or input values) using the formula:

weight_final = weight_initial / (scale * sky)

where the sky values are those from the image statistics before conversion to zero level shifts and adjustment to zero mean over all images. The reasoning is that if the zero level is high the sky brightness is high and so the S/N is lower and the weight should be lower. If any sky value determined from the image statistics comes out to be negative a warning is given and the none of the weight are adjusted for sky levels.

The weights are not adjusted when the zero offsets are input from a file or keyword since these values do not imply the actual image sky value. In this case if one wants to account for different sky statistics in the weights the user must specify the weights in a file taking explicit account of changes in the weights due to different sky statistics.

PIXEL MASKS

A pixel mask is a type of IRAF file having the extension ".pl" which identifies an integer value with each pixel of the images to which it is applied. The integer values may denote regions, a weight, a good or bad flag, or some other type of integer or integer bit flag. In the common case where many values are the same this file is compacted to be small and efficient to use. It is also most compact and efficient if the majority of the pixels have a zero mask value so frequently zero is the value for good pixels. Note that these files, while not stored as a strict pixel array, may be treated as images in programs. This means they may be created by programs such as mkpattern, edited by imedit, examined by imexamine, operated upon by imarith, graphed by implot, and displayed by display.

At the time of introducing this task, generic tools for creating pixel masks have yet to be written. There are two ways to create a mask in V2.10. First if a regular integer image can be created then it can be converted to pixel list format with imcopy:

cl> imcopy template plfile.pl

by specifically using the .pl extension on output. Other programs that can create integer images (such mkpattern or ccdred.badpiximage) can create the pixel list file directly by simply using the ".pl" extension in the output image name.

To use pixel masks with combine one must associate a pixel mask file with an image by entering the pixel list file name in the image header under the keyword BPM (bad pixel mask). This can be done with hedit. Note that the same pixel mask may be associated with more than one image as might be the case if the mask represents defects in the detector used to obtain the images.

If a pixel mask is associated with an image the mask is used when the masktype parameter is set to a value other than "none". Note that when it is set to "none" mask information is not used even if it exists for the image. The values of masktype which apply masks are "goodvalue", "badvalue", "goodbits", and "badbits". They are used in conjunction with the maskvalue parameter. When the mask type is "goodvalue" the pixels with mask values matching the specified value are included in combining and all others are rejected. Similarly, for a mask type of "badvalue" the pixels with mask values matching the specified value are rejected and all others are accepted. The bit types are useful for selecting a combination of attributes in a mask consisting of bit flags. The mask value is still an integer but is interpreted by bitwise comparison with the values in the mask file.

If a mask operation is specified and an image has no mask image associated with it then the mask values are taken as all zeros. In those cases be careful that zero is an accepted value otherwise the entire image will be rejected.

In the case of combining the higher dimensions of an image into a lower dimensional image, the "project" option, the same pixel mask is applied to all of the data being combined; i.e. the same 2D pixel mask is applied to every plane of a 3D image. This is because a higher dimensional image is treated as a collection of lower dimensional images having the same header and hence the same bad pixel mask. It would be tempting to use a bad pixel mask with the same dimension as the image being projected but this is not currently how the task works.

When the number of input images exceeds the maximum number of open files allowed by IRAF (currently about 115 .imh or 57 .hhh images) the input images are stacked and combined with the project option. This means that the bad pixel mask from the first input image will be applied to all the images.

THRESHOLD REJECTION

In addition to rejecting masked pixels, pixels in the unscaled input images which are below or above the thresholds given by the parameters lthreshold and hthreshold are rejected. Values of INDEF mean that no threshold value is applied. Threshold rejection may be used to exclude very bad pixel values or as an alternative way of masking images. In the latter case one can use a task like imedit or imreplace to set parts of the images to be excluded to some very low or high magic value.

REJECTION ALGORITHMS

The reject parameter selects a type of rejection operation to be applied to pixels not masked or thresholded. If no rejection operation is desired the value "none" is specified.

MINMAX A specified fraction of the highest and lowest pixels are rejected. The fraction is specified as the number of high and low pixels, the nhigh and nlow parameters, when data from all the input images are used. If pixels have been rejected by offseting, masking, or thresholding then a matching fraction of the remaining pixels, truncated to an integer, are used. Thus,

nl = n * nlow/nimages + 0.001
nh = n * nhigh/nimages + 0.001

where n is the number of pixels surviving offseting, masking, and thresholding, nimages is the number of input images, nlow and nhigh are task parameters and nl and nh are the final number of low and high pixels rejected by the algorithm. The factor of 0.001 is to adjust for rounding of the ratio.

As an example with 10 input images and specifying one low and two high pixels to be rejected the fractions to be rejected are nlow=0.1 and nhigh=0.2 and the number rejected as a function of n is:

n   0  1  2  3  4  5  6  7  8  9 10
nl  0  0  0  0  0  0  0  0  0  0  1
nh  0  0  0  0  0  1  1  1  1  1  2

CCDCLIP If the images are obtained using a CCD with known read out noise, gain, and sensitivity noise parameters and they have been processed to preserve the relation between data values and photons or electrons then the noise characteristics of the images are well defined. In this model the sigma in data values at a pixel with true value <I>, as approximated by the median or average with the lowest and highest value excluded, is given by:

sigma = ((rn / g) ** 2 + <I> / g + (s * <I>) ** 2) ** 1/2

where rn is the read out noise in electrons, g is the gain in electrons per data value, s is a sensitivity noise given as a fraction, and ** is the exponentiation operator. Often the sensitivity noise, due to uncertainties in the pixel sensitivities (for example from the flat field), is not known in which case a value of zero can be used. See the task stsdas.wfpc.noisemodel for a way to determine these vaues (though that task expresses the read out noise in data numbers and the sensitivity noise parameter as a percentage).

The read out noise is specified by the rdnoise parameter. The value may be a numeric value to be applied to all the input images or a image header keyword containing the value for each image. Similarly, the parameter gain specifies the gain as either a value or image header keyword and the parameter snoise specifies the sensitivity noise parameter as either a value or image header keyword.

The algorithm operates on each output pixel independently. It starts by taking the median or unweighted average (excluding the minimum and maximum) of the unrejected pixels provided there are at least two input pixels. The expected sigma is computed from the CCD noise parameters and pixels more that lsigma times this sigma below or hsigma times this sigma above the median or average are rejected. The process is then iterated until no further pixels are rejected. If the average is used as the estimator of the true value then after the first round of rejections the highest and lowest values are no longer excluded. Note that it is possible to reject all pixels if the average is used and is sufficiently skewed by bad pixels such as cosmic rays.

If there are different CCD noise parameters for the input images (as might occur using the image header keyword specification) then the sigmas are computed for each pixel from each image using the same estimated true value.

If the images are scaled and shifted and the sigscale threshold is exceedd then a sigma is computed for each pixel based on the image scale parameters; i.e. the median or average is scaled to that of the original image before computing the sigma and residuals.

After rejection the number of retained pixels is checked against the nkeep parameter. If there are fewer pixels retained than specified by this parameter the pixels with the smallest residuals in absolute value are added back. If there is more than one pixel with the same absolute residual (for example the two pixels about an average or median of two will have the same residuals) they are all added back even if this means more than nkeep pixels are retained. Note that the nkeep parameter only applies to the pixels used by the clipping rejection algorithm and does not apply to threshold or bad pixel mask rejection.

This is the best clipping algorithm to use if the CCD noise parameters are adequately known. The parameters affecting this algorithm are reject to select this algorithm, mclip to select the median or average for the center of the clipping, nkeep to limit the number of pixels rejected, the CCD noise parameters rdnoise, gain and snoise, lsigma and hsigma to select the clipping thresholds, and sigscale to set the threshold for making corrections to the sigma calculation for different image scale factors.

CRREJECT This algorithm is identical to "ccdclip" except that only pixels above the average are rejected based on the hsigma parameter. This is appropriate for rejecting cosmic ray events and works even with two images.

SIGCLIP The sigma clipping algorithm computes at each output pixel the median or average excluding the high and low values and the sigma about this estimate. There must be at least three input pixels, though for this method to work well there should be at least 10 pixels. Values deviating by more than the specified sigma threshold factors are rejected. These steps are repeated, except that after the first time the average includes all values, until no further pixels are rejected or there are fewer than three pixels.

After rejection the number of retained pixels is checked against the nkeep parameter. If there are fewer pixels retained than specified by this parameter the pixels with the smallest residuals in absolute value are added back. If there is more than one pixel with the same absolute residual (for example the two pixels about an average or median of two will have the same residuals) they are all added back even if this means more than nkeep pixels are retained. Note that the nkeep parameter only applies to the pixels used by the clipping rejection algorithm and does not apply to threshold or bad pixel mask rejection.

The parameters affecting this algorithm are reject to select this algorithm, mclip to select the median or average for the center of the clipping, nkeep to limit the number of pixels rejected, lsigma and hsigma to select the clipping thresholds, and sigscale to set the threshold for making corrections to the sigma calculation for different image scale factors.

AVSIGCLIP The averaged sigma clipping algorithm assumes that the sigma about the median or mean (average excluding the low and high values) is proportional to the square root of the median or mean at each point. This is described by the equation:

sigma(column,line) = sqrt (gain(line) * signal(column,line))

where the estimated signal is the mean or median (hopefully excluding any bad pixels) and the gain is the estimated proportionality constant having units of photons/data number.

This noise model is valid for images whose values are proportional to the number of photons recorded. In effect this algorithm estimates a detector gain for each line with no read out noise component when information about the detector noise parameters are not known or available. The gain proportionality factor is computed independently for each output line by averaging the square of the residuals (at points having three or more input values) scaled by the median or mean. In theory the proportionality should be the same for all rows but because of the estimating process will vary somewhat.

Once the proportionality factor is determined, deviant pixels exceeding the specified thresholds are rejected at each point by estimating the sigma from the median or mean. If any values are rejected the median or mean (this time not excluding the extreme values) is recomputed and further values rejected. This is repeated until there are no further pixels rejected or the number of remaining input values falls below three. Note that the proportionality factor is not recomputed after rejections.

If the images are scaled differently and the sigma scaling correction threshold is exceedd then a correction is made in the sigma calculations for these differences, again under the assumption that the noise in an image scales as the square root of the mean intensity.

After rejection the number of retained pixels is checked against the nkeep parameter. If there are fewer pixels retained than specified by this parameter the pixels with the smallest residuals in absolute value are added back. If there is more than one pixel with the same absolute residual (for example the two pixels about an average or median of two will have the same residuals) they are all added back even if this means more than nkeep pixels are retained. Note that the nkeep parameter only applies to the pixels used by the clipping rejection algorithm and does not apply to threshold or bad pixel mask rejection.

This algorithm works well for even a few input images. It works better if the median is used though this is slower than using the average. Note that if the images have a known read out noise and gain (the proportionality factor above) then the "ccdclip" algorithm is superior. The two algorithms are related in that the average sigma proportionality factor is an estimate of the gain.

The parameters affecting this algorithm are reject to select this algorithm, mclip to select the median or average for the center of the clipping, nkeep to limit the number of pixels rejected, lsigma and hsigma to select the clipping thresholds, and sigscale to set the threshold for making corrections to the sigma calculation for different image scale factors.

PCLIP The percentile clipping algorithm is similar to sigma clipping using the median as the center of the distribution except that, instead of computing the sigma of the pixels from the CCD noise parameters or from the data values, the width of the distribution is characterized by the difference between the median value and a specified "percentile" pixel value. This width is then multipled by the scale factors lsigma and hsigma to define the clipping thresholds above and below the median. The clipping is not iterated.

The pixel values at each output point are ordered in magnitude and the median is determined. In the case of an even number of pixels the average of the two middle values is used as the median value and the lower or upper of the two is the median pixel when counting from the median pixel to selecting the percentile pixel. The parameter pclip selects the percentile pixel as the number (if the absolute value is greater than unity) or fraction of the pixels from the median in the ordered set. The direction of the percentile pixel from the median is set by the sign of the pclip parameter with a negative value signifying pixels with values less than the median. Fractional values are internally converted to the appropriate number of pixels for the number of input images. A minimum of one pixel and a maximum corresponding to the extreme pixels from the median are enforced. The value used is reported in the log output. Note that the same percentile pixel is used even if pixels have been rejected by offseting, masking, or thresholding; for example, if the 3nd pixel below the median is specified then the 3rd pixel will be used whether there are 10 pixels or 5 pixels remaining after the preliminary steps.

Some examples help clarify the definition of the percentile pixel. In the examples assume 10 pixels. The median is then the average of the 5th and 6th pixels. A pclip value of 2 selects the 2nd pixel above the median (6th) pixel which is the 8th pixel. A pclip value of -0.5 selects the point halfway between the median and the lowest pixel. In this case there are 4 pixels below the median, half of that is 2 pixels which makes the percentile pixel the 3rd pixel.

The percentile clipping algorithm is most useful for clipping small excursions, such as the wings of bright objects when combining disregistered observations for a sky flat field, that are missed when using the pixel values to compute a sigma. It is not as powerful, however, as using the CCD noise parameters (provided they are accurately known) to clip about the median.

The parameters affecting this algorithm are reject to select this algorithm, pclip to select the percentile pixel, nkeep to limit the number of pixels rejected, and lsigma and hsigma to select the clipping thresholds.

GROW REJECTION

Neighbors of pixels rejected by the rejection algorithms along image lines may also be rejected. The number of neighbors to be rejected on either side is specified by the grow parameter. The rejection only applies to neighbors along each image line. This is because the task operates independently on each image line and does not have the ability to go back to previous lines or maintain a list of rejected pixels to later lines.

This rejection step is also checked against the nkeep parameter and only as many pixels as would not violate this parameter are rejected. Unlike it's application in the rejection algorithms at this stage there is no checking on the magnitude of the residuals and the pixels retained which would otherwise be rejected are randomly selected.

COMBINING

After all the steps of offsetting the input images, masking pixels, threshold rejection, scaling, and applying a rejection algorithms the remaining pixels are combined and output. The pixels may be combined by computing the median or by computing a weighted average.

SIGMA OUTPUT

In addition to the combined image and optional sigma image may be produced. The sigma computed is the standard deviation, corrected for a finite population by a factor of n/(n-1), of the unrejected input pixel values about the output combined pixel values.

Examples

1. To average and median images without any other features:

cl> combine obj* avg combine=average reject=none
cl> combine obj* med combine=median reject=none

2. To reject cosmic rays:

cl> combine obs1,obs2 Obs reject=crreject rdnoise=5.1, gain=4.3

3. To make a grid for display purposes with 21 64x64 images:

cl> combine @list grid offset="grid 5 65 5 65"

4. To apply a mask image with good pixels marked with a zero value and bad pixels marked with a value of one:

cl> hedit ims* bpm badpix.pl add+ ver-
cl> combine ims* final combine=median masktype=goodval

5. To scale image by the exposure time and then adjust for varying sky brightness and make a weighted average:

cl> combine obj* avsig combine=average reject=avsig \
>>> scale=exp zero=mode weight=exp  expname=exptime

Time requirements

The following times were obtain with a Sun 4/470. The tests combine 1000x200 images consisting of Poisson noise and cosmic rays generated with the artdata package. The times, especially the total time, are approximate and depend on user loads.

IMAGES:   Number of images (1000x200) and datatype (R=real, S=short)
COMBINE:  Combine option
REJECT:   Rejection option with grow = 0
              minmax:    nlow = 1, nhigh = 1
              ccdclip:   lsigma = 3., hsigma = 3, sigscale = 0.
              sigclip:   lsigma = 3., hsigma = 3, sigscale = 0.
              avsigclip: lsigma = 3., hsigma = 3, sigscale = 0.
              pclip:     lsigma = 3., hsigma = 3, pclip = -0.5
              /a:        mclip = no  (clip about the average)
              /m:        mclip = yes (clip about the median)
O M T S:  Features used (Y=yes, N=no)
O:        offset = "grid 5 10 2 10"
M:        masktype = goodval, maskval = 0
              Pixel mask has 2 bad lines and 20 bad columns
T:        lthreshold = INDEF, hthreshold = 1100.
S:        scale = mode, zero = none, weight = mode
TIME:     cpu time in seconds, total time in minutes and seconds

IMAGES  COMBINE  REJECT        O M T S     TIME

  10R   average  none          N N N N    1.3 0:08
  10R   average  minmax        N N N N    4.3 0:10
  10R   average  pclip         N N N N   17.9 0:32
  10R   average  ccdclip/a     N N N N   11.6 0:21
  10R   average  crreject/a    N N N N   11.4 0:21
  10R   average  sigclip/a     N N N N   13.6 0:29
  10R   average  avsigclip/a   N N N N   15.9 0:35
  10R   average  ccdclip/m     N N N N   16.9 0:32
  10R   average  crreject/m    N N N N   17.0 0:28
  10R   average  sigclip/m     N N N N   19.6 0:42
  10R   average  avsigclip/m   N N N N   20.6 0:43

  10R   median   none          N N N N    6.8 0:17
  10R   median   minmax        N N N N    7.8 0:15
  10R   median   pclip         N N N N   16.9 1:00
  10R   median   ccdclip/a     N N N N   18.0 0:34
  10R   median   crreject/a    N N N N   17.7 0:30
  10R   median   sigclip/a     N N N N   21.1 1:13
  10R   median   avsigclip/a   N N N N   23.1 0:41
  10R   median   ccdclip/m     N N N N   16.1 0:27
  10R   median   crreject/m    N N N N   16.0 0:27
  10R   median   sigclip/m     N N N N   18.1 0:29
  10R   median   avsigclip/m   N N N N   19.6 0:32

  10R   average  none          N N N Y    6.1 0:36
  10R   median   none          N N N Y   10.4 0:49
  10R   median   pclip         N N N Y   20.4 1:10
  10R   median   ccdclip/m     N N N Y   19.5 0:36
  10R   median   avsigclip/m   N N N Y   23.0 1:06

  10R   average  none          N Y N N    3.5 0:12
  10R   median   none          N Y N N    8.9 0:21
  10R   median   pclip         N Y N N   19.9 0:45
  10R   median   ccdclip/m     N Y N N   18.0 0:44
  10R   median   avsigclip/m   N Y N N   20.9 0:28

  10R   average  none          Y N N N    4.3 0:13
  10R   median   none          Y N N N    9.6 0:21
  10R   median   pclip         Y N N N   21.8 0:54
  10R   median   ccdclip/m     Y N N N   19.3 0:44
  10R   median   avsigclip/m   Y N N N   22.8 0:51

  10R   average  none          Y Y Y Y   10.8 0:22
  10R   median   none          Y Y Y Y   16.1 0:28
  10R   median   pclip         Y Y Y Y   27.4 0:42
  10R   median   ccdclip/m     Y Y Y Y   25.5 0:39
  10R   median   avsigclip/m   Y Y Y Y   28.9 0:44

  10S   average  none          N N N N    2.2 0:06
  10S   average  minmax        N N N N    4.6 0:12
  10S   average  pclip         N N N N   18.1 0:33

Revisions

COMBINE V2.11: The limit of the number of images that may be combined has been removed. If the number of images exceeds the maximum number of open images permitted then the images are stacked in a single temporary image and then combined with the project option. Note that this will double the amount of diskspace temporarily. There is also a limitation in this case that the bad pixel mask from the first image in the list will be applied to all the images. Integer offsets may be determined from the image world coordinate system.

COMBINE V2.10.3: The output pixel datatype parameter, outtype was previously ignored and the package pixeltype was used. The task output pixel type parameter is now used. The factors specified by an @file or keyword are not normalized.

COMBINE V2.10.2: The weighting was changed from using the square root of the exposure time or image statistics to using the values directly. This corresponds to variance weighting. Other options for specifying the scaling and weighting factors were added; namely from a file or from a different image header keyword. The nkeep parameter was added to allow controlling the maximum number of pixels to be rejected by the clipping algorithms. The snoise parameter was added to include a sensitivity or scale noise component to the noise model. Errors will now delete the output images.

COMBINE V2.10

This task was greatly revised to provide many new features. These features are:

o Bad pixel masks
o Combining offset and different size images
o Blank value for missing data
o Combining across the highest dimension (the project option)
o Separating threshold rejection, the rejection algorithms,
  and the final combining statistic
o New CCDCLIP, CRREJECT, and PCLIP algorithms
o Rejection now may reject more than one pixel per output pixel
o Choice of a central median or average for clipping
o Choice of final combining operation
o Simultaneous multiplicative and zero point scaling

Limitations

Though the previous limit on the number of images that can be combined was removed in V2.11 the method has the limitation that only a single bad pixel mask will be used for all images.