path: root/pkg/images/immatch/doc/imcombine.hlp

.help imcombine Aug01 images.immatch
.ih
NAME
imcombine -- Combine images using various algorithms
.ih
USAGE	
imcombine input output
.ih
PARAMETERS
.ls input
List of input images to combine.  If the \fIproject\fR parameter is "no"
then all input images must have the same dimensionality though they may
be of different sizes.  Otherwise each input image is handled separately
and they may have different dimensionalities.
.le


When the \fIproject\fR parameter is "no" all the input images are combined
into a single output file.  In this case the following parameters specify
only a single file name.  Otherwise each input image is combined by
projecting (combining across) the highest dimension to produce a lower
dimensional image.  For this type of combining there is one output for each
input and so the following parameters specify matching lists.

.ls output
Output combined image(s).  If there are fewer than 100 input images the
names of the input images are recorded in header keywords IMCMBnnn.
.le
.ls headers = "" (optional)
Optional output multiextension FITS file(s).  The extensions are dataless
headers from each input image.
.le
.ls bpmasks = "" (optional)
Optional output bad pixel mask(s) with good values of 0 and bad values of
1.  Output pixels are marked as bad when no input pixels contributed to the
output pixel.  The file name is also added to the output image header under
the keyword BPM.
.le
.ls rejmask = "" (optional)
Optional output mask file(s) identifying rejected or excluded pixels.  The
pixel mask is the size of the output image but there is one extra dimension
with length equal to the number of input images.  Each element of the
highest dimension is a mask corresponding to an input image with values of
1 for rejected or excluded pixels and values of 0 for pixels which were
used.  The order of the masks is the order of the input images and image
header keywords, indexed by the pixel coordinate of the highest dimension
identify the input images.  Note that the pixel positions are in the output
pixel coordinate system.
.le
.ls nrejmasks = "" (optional)
Optional output pixel mask(s) giving the number of input pixels rejected or
excluded from the input images.
.le
.ls expmasks = "" (optional)
Optional output exposure mask(s) giving the sum of the exposure values of
the input images with non-zero weights that contributed to that pixel.
Since masks are integer, the exposure values may be scaled to preserve
dynamic range and fractional significance.  The scaling values are given in
the header under the keywords MASKSCAL and MASKZERO.  Exposure values are
computed from the mask values by scale * value + zero where scale is the
value of the MASKSCAL keyword and zero is the value of the MASKZERO
keyword.
.le
.ls sigma = "" (optional)
Optional output sigma image(s).  The sigma is the standard deviation,
corrected for a finite population, of the input pixel values (excluding
rejected pixels) about the output combined pixel values.
.le

.ls imcmb = "$I" (optional)
A keyword in the input images that is copied
to one of the IMCMBnnn keywords in the output image.  A null string
does not set the IMCMBnnn keywords nor deletes any existing keywords.
Any other value will delete existing keywords before creating new ones.
The special value "$I" specifies the basename of the input image name.
If a keyword is specified that does not exist in the input image(s) then
no ICMB keyword will be produced; it is not a error for the keyword to
not exist.
.le
.ls logfile = "STDOUT" (optional)
Optional output log file.  If no file is specified then no log information is
produced.  The special filename "STDOUT" prints log information to the
terminal.
.le

.ls combine = "average" (average|median|lmedian|sum|quadrature|nmodel)
Type of combining operation performed on the final set of pixels (after
offsetting, masking, thresholding, and rejection).  The choices are:

.nf
    average - weighted average
     median - median
    lmedian - median except use the lower value if only two
        sum - (weighted) sum
 quadrature - weighted quadrature average
     nmodel - weighted quadrature average of noise model values
.fi

The details of each choice is given in the DESCRIPTION.
Note that if weights are used then the weighted "sum" is the same as
the weighted "average" since the weights are normalized to unit total weight.
The "lmedian" option is intended for minimizing the effects of cosmic rays
when there are more than two images but some pixels may only have two
contributing images.  The "quadrature" and "nmodel" options are used
for error propagation either with input sigma images (quadrature) or where the
pixel sigmas may be computed by the noise model used by this task (nmodel).
.le
.ls reject = "none" (none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip)
Type of rejection operation performed on the pixels remaining after offsetting,
masking and thresholding.  The algorithms are described in the
DESCRIPTION section.  The rejection choices are:

.nf
      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
.fi

.le
.ls 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.
.le
.ls outtype = "real" (none|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.  When there is a mixture of
short and unsigned short images the highest precedence become integer.
The datatypes may be abbreviated to a single character.
.le
.ls outlimits = ""
Output region limits specified as pairs of whitespace separated values.
The first two numbers are the limits along the first output image dimension,
the next two numbers are the limits along the second dimension, and so on.
If the higher dimension limits are not specified they default to the full
range.  Therefore, if no limits are specified then the full output is
created.  Note that the output size is computed from all the input images
including offsets if specified and the coordinates are relative to that
size.
.le
.ls offsets = "none" (none|wcs|world|physical|grid|<filename>)
Integer offsets to add to each image axes.  The options are:
.ls "none"
No offsets are applied.
.le
.ls "wcs" or "world"
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.
.le
.ls "physical"
The physical coordinate system defined by the IRAF LTM/LTV keywords
define the offsets.
.le
.ls "grid"
A uniform grid of offsets is specified by a string of the form

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

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.
.le
.ls <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.
.le
.le
.ls masktype = "none"
Type of pixel masking to use.  The choices are

.nf
         none - No pixel masking
    goodvalue - good pixels defined by maskvalue parameter
     badvalue - bad pixels defined by maskvalue parameter
      novalue - pixels with no value defined by maskvalue parameter
     goodbits - good pixels defined by maskvalue parameter
      badbits - bad pixels defined by maskvalue parameter
.fi

Except for "none", these choices use the mask specified by the header
keyword BPM.  To use a different keyword to specify the mask the syntax
is

.nf
  !<keyword> [goodvalue|badvalue|novalue|goodbits|badbits]
.fi

where if the optional second word is missing the default is "goodvalue".

If "none" (or "") no pixel masking is done
even if an image has an associated  pixel mask.  The masking defines
pixels to be used (good) and not used (bad).  The types use the
"maskvalue" parameter to define a single value (either as a number or
set of bits) for good or bad and all other values are treated as the
opposite; i.e. bad or good respectively.

The "novalue" choice uses 0 as the good value and all other values are
bad.  However, the "maskvalue" parameter defines a mask value for pixels
with no value such as occurs from rebinning at the edges or stacking where
there is no overlap at all.  The distinction pixels is that when a final pixel
has no overlapping data because all input pixels have a "no value" flag
the blank value is output while if there is no good data then pixels which
are have other than the "no value" flag are combined as if they were good
to produce a representative output value.  An output mask will have a
value of 0 for pixels where at least one good input value was present,
a value of 1 when there was no overlapping data, and a value of 2 when
bad data was used.
.le
.ls maskvalue = 0
Mask value used with the \fImasktype\fR 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.
.le
.ls blank = 0.
Output value to be used when there are no pixels for combining after any
rejection.
.le

.ls scale = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>)
Multiplicative image scaling to be applied.  The choices are none, multiply
by the reciprocal of the mode, median, or mean of the specified statistics
section, multiply by the reciprocal of the exposure time in the image header,
multiply by the values in a specified file, or multiply 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.
.le
.ls zero = "none" (none|mode|median|mean|@<file>|!<keyword>)
Additive zero level image shifts to be applied.  The choices are none, add
the negative of the mode, median, or mean of the specified statistics
section, add the values given in a file, or add the 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.
.le
.ls 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.
.le
.ls 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.
.le
.ls  expname = ""
Image header keyword to be used with the exposure scaling and weighting
options.  Also if an exposure keyword is specified that keyword will be
added to the output image using a weighted average of the input exposure
values.
.le

.ce
Algorithm Parameters
.ls 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.
.le
.ls 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 non-overlap, then the fraction of the remaining pixels, truncated
to an integer, is used.
.le
.ls 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.  The latter is in addition to pixels
missing due to non-overlapping offsets, bad pixel masks, or thresholds.
.le
.ls 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.
.le
.ls 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.
.le
.ls rdnoise = "0.", gain = "1.", snoise = "0." (ccdclip, crreject)
Readout noise in electrons, gain in electrons/DN, and sensitivity noise as
a fraction.  These parameters are used with the "ccdclip" and "crreject"
algorithms as well as with the "nmodel" combining option.  The values may
be either numeric or an image header keyword which contains the value.
The noise model for a pixel is:

.nf
    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
.fi

where DN is the data number and Ne is the number of electrons.  Sensitivity
noise typically comes from noise introduced during flat fielding.
.le
.ls 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.
.le
.ls 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.
.le
.ls grow = 0.
Radius in pixels for additional pixel to be rejected in an image with a
rejected pixel from one of the rejection algorithms.  This applies only to
pixels rejected by one of the rejection algorithms and not the masked or
threshold rejected pixels.
.le

.ce
Environment Variables

.ls imcombine_maxmemory (default = 250000000)
This task tries to use the maximum possible memory for efficiency when
dealing with lots of data and is designed to reduce memory usage if
memory allocation fails.  However, there may be cases where this adjustment
fails so this variable allows forcing the task to stay within a smaller
allocation.  This variable is in bytes and the default is the amount
generally returned by the system.  It is large because of virtual memory
functionality.  If problems are encountered one should try setting this
variable to a smaller size until, hopefully, the out of memory errors
disappear.
.le
.ls imcombine_option (default = 1)
This environment variable is used to select certain experimental or
diagnostic options.  If this variable has the value 1, the default when the
variable is undefined, then when the number of images exceeds the number of
files that can be kept open under IRAF (currently this means more than 4000
images) the images are closed and opened as needed.  This is in contrast to
the previous method, when the variable has the value 0, which first builds
a single stacked image of a higher dimension from the input images.  This
method requires the images all have the same size and also that there be
sufficient disk space for the stacked image and that the image  be less
than 2Gb in size.
.le
.ih
DESCRIPTION
A set of images or the highest dimension elements (for example
the planes in an image cube) are combined by weighted averaging, medianing,
or summing.  Pixels may be rejected from the combining by using pixel
masks, threshold levels, and rejection algorithms.  The images may be
scaled, before rejections, multiplicatively, additively, or both based on
image statistics, image header keywords, or text files.  The images may be
combined with integer pixel coordinate offsets, possibly determined using
the world coordinate system of the images, to produce an image bigger than
any of the input images.

The input images to be combined are specified by a list.  If the
\fBproject\fR 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.  If \fBproject\fR is "no" then the entire input list is combined to
form a single output image.   In this case the images must all have the
same dimensionality but they may have different sizes.  There is a software
limit of approximately 4000 images which may be open
simultaneously.  To combine more than this number the program may either
create a temporary stacked image, requiring the images to be of the same
size, or repeatedly open and close the images.  See the "Environment
Variables" in the PARAMETERS section.

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 value of a keyword in the input images, where the
keyword is specified by the parameter \fIimcmb\fR, is written to an
indexed keyword IMCMBnnn.  The purpose of the ICMBnnn keywords is to
identify the contributors to the output image.  One common choice is
the input image name though other identifiers may be used.

If a bad pixel mask is created, the name of the mask will be included in the
output image header under the keyword BPM.  The output pixel type is set by
the parameter \fIouttype\fR.  If left blank then the input datatype of
highest precision is used.  If there is a mixture of short and unsigned
short images then the highest precision is integer.

In addition to one or more output combined images there are some optional
output files which may be specified as described in the OPTIONAL OUTPUT
section.

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

.nf
o   Check the input images and stack them if needed
o   Set the input image offsets and the final output image size.
o   Set the input image scales and weights possibly by computing
    image statistics
o   Write the log file and optional header output
.fi

For each output image line:

.nf
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 and other optional images.
.fi

OPTIONAL OUTPUTS

There are various additional outputs that may be produced by providing
the filenames.

.ls Headers
The output image can only have one set of header keywords which are
inherited from the first input image in the input list.  Copies of all the
input headers may be stored in a multiextension FITS file specified by the
\fIheaders\fR parameter.  The extension names are the input image names.
The extensions are dataless headers.  Since this means the image sizes are
lost, AXLEN keywords are added.  Also the keywords INIMAGE and OUTIMAGE are
added giving the name of the input image and the name of the output
combined image.
.le
.ls Bad Pixel Masks
The \fIbpmasks\fR parameter produces optional output bad pixel mask(s) with
good values of 0 and bad values of 1.  Output pixels are marked as bad when
no input pixels contributed to the output pixel.  The file name is also
added to the output image header under the keyword BPM.
.le
.ls Rejection Masks
The \fIrejmasks\fR parameter produces optional output mask file(s)
identifying rejected or excluded pixels.  The pixel mask is the size of the
output image.  There is one extra dimension with length equal to the number
of input images.  Each element of the highest dimension is a mask
corresponding to an input image with values of 1 for rejected or excluded
pixels and values of 0 for pixels which were used.  The order of the masks
is the order of the input images and image header keywords indexed by the
element identify the input images.  Note that the pixel positions are in
the output pixel coordinate system.

This mask is the only way to record whether a particular input image pixel
contributed to the output image.  As an example, consider the case of
three input two dimensional images of sizes 1020x1020, 1010x1010, and
1000x1000 with relative offsets of (0,0), (10,10), and (20,20).  The output
image would then be 1020x1020.

Suppose that the only input pixels not used are pixels (1,1) in each input
image.  Because of the offsets the first 10 rows and columns of the output
will be based on just a single pixel except for (1,1) which has no input
pixels.  The next 10 rows and columns of the output will be a combination
of 2 input pixels except (11,11) which is just based on pixel (11,11)
in the first input image.  Finally all other pixels except (21,21) will be
a combination of 3 values.

The rejection mask will be three dimensional of size 1020x1020x3.  Plane 1
will correspond to the first input image, plane 2 to the second, and so
on.  All of the pixels will be zero except for the following pixels
which will have a value of 1. In the first plane only pixel (1,1,1) will be
one.  In the second plane the first 10 rows and columns and pixel (11,11,2)
will be one.  And in the third plane, the first 20 rows and columns and pixel
(21,21,3) will be one.  So if we want to know about output pixel (11,11)
the rejection mask will tell us that pixels from the second and third
images were excluded.

This is a complex example because of the offsets and dissimilar sizes.
In the more common and simpler case of equal sizes and registered images,
the mask
planes would have one to indicate that the pixel in the input image at
that coordinate was not used.  For instance if pixel (12,15,2) is one
in the rejection mask then pixel (12,15) in the second input image was
excluded.

Note that one can use image sections to extract a mask matching the input
image.  For the example case with the offsets masks in the input
coordinates can be extracted with the commands

.nf
    cl> imcopy rejmask[*,*,1] mask1
    cl> imcopy rejmask[11:1020,11:1020,2] mask2
    cl> imcopy rejmask[21:1020,21:1020,3] mask3
.fi

For the case of equal sized and registered images one could also use
\fBimslice\fR.
.le
.ls Mask of the Number of Rejected Pixels
The \fInrejmasks\fR parameter produces optional pixel mask(s) giving the
number of input pixels rejected or excluded from the input images.  This is
equivalent to projecting the rejection mask described previously by summing
along the highest dimension.  Note that in this mask a value of 0 indicates
all the input pixels were used to create the output pixel and a value equal
to the number of input images indicate no input pixels were used.
.le
.ls Exposure Masks
The \fIexpmasks\fR parameter produces optional output exposure mask(s)
giving the sum of the exposure values of the input images with non-zero
weights that contributed to that pixel.  Since masks are integer, the
exposure values may be scaled to preserve dynamic range and fractional
significance.  The scaling values are given in the header under the
keywords MASKSCAL and MASKZERO.  Exposure values are computed from the mask
values by scale * value + zero where scale is the value of the MASKSCAL
keyword and zero is the value of the MASKZERO keyword.
.le
.ls Sigma of Combined Pixels
The \fIsigma\fR parameter produces optional output sigma image(s).  The
sigma is the standard deviation, corrected for a finite population, of the
input pixel values (excluding rejected pixels) about the output combined
pixel values.
.le
.ls Output Log File
The \fIlogfile\fR parameter produces an optional output log file.  If no
file is specified then no log information is produced.  The special
filename "STDOUT" prints log information to the terminal.
.le

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 \fIproject\fR 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 columns, the first two images in the next 10 columns and
all three images starting in the 21st column.  At the 21st output column
the 21st column of the first image will be combined with the 11th 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 120
columns corresponding to the 20 column offset in the third image.
Note that this same output image size is computed and used as the
basis for the \fIoutlimits\fR parameter to specify a subregion to
actually be output.

The input image offsets are set using the \fIoffset\fR 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

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

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
non-overlapping 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 \fIblank\fR parameter between the images.

The offsets may be defined in a file by specifying the file name
in the \fIoffset\fR 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:

.nf
	0 0
	10 0
	20 0
.fi

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:

.nf
	225 15
	235 15
	245 15
.fi

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

.nf
	# Absolute
	10 10
	20 10
	30 10
.fi

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.

OUTPUT OF SUBREGIONS

The output image size is computed from all of the input images including
the offsets as described previously.  The \fIoutlimits\fR may be used to
specify a subregion of this full size to be created.  The syntax of this
parameter is pairs of whitespace separated numbers selecting the first and last
pixel in each output dimension.  The pairs for each dimension are also
whitespace separated and are given in order of the dimensions.  Any missing
values at the end of the string default to the full limits of the output
image.  If the requested limits fall outside the full output image they are
reset to the size of the full computed output size.

As an example, consider combining 10 images of size 1000x1000 with offsets
of 0, 1, ..., 9 along the first dimension.  Because of the offsets the full
output size is 1010x1000.  To output only the region [1:100,101:200]
of this full size the parameter value would be the string "1 100 101 200".
Note that if the value was just "1 100" then the output region would
be [1:100,1:1000].

The intended purpose for this option is to allow creating subregions using
a smaller number of images in the case of offset data taken at a raster of
positions.  This is important since when the number of images becomes too
large (>4000) the program either has to prestack the images into a higher
dimensional single image (requiring equal sized images) or utilize an
inefficient algorithm where images are opened and closed for each input
line.  A detail of how this task works is that it is the number of images
required to be accessed for each output line that is significant.

The following example was developed when the maximum number of images
open at one time was ~240.  In V2.12 the number was increased to
more than 4000 so it is not as relevant though it may apply to very
large surveys with many small images.

As an example, consider a survey of a region of the sky composed of 8000
images which are each 500x1000.  The offsets between each image are 50
pixels along the first dimension and 900 pixels along the second dimension,
give or take a few pixels due to telescope pointing errors.  Thus this
survey consists of strips of exposures.  Within a strip the images over by
about 450 pixels.  Between strips the overlap is 100 pixels.  So the
strips consist 400 exposures and there are 20 strips.

The full size of this survey is then about 20450x18900.  At any point in a
single strip the number of images contributing is no more than 10.
Including the overlap of the strips the maximum number is then 20.  In
order to combine the data for such a survey one would like to create
subregion outputs which are 120 images from each strip.  The lines where
the two strips overlap then require 240 images.  To produce roughly equal
size regions we choose sizes along the first dimension of 5200 pixels.  The
number of lines in the output subregions might be the full size of the
survey.  However, it might be desirable to also  break the second dimension
into blocks for ease of display and manipulation.

The method for combining this example survey is then to combine the data in
four groups along the first dimension to produce subimages each 5200 pixels
wide which have no overlap.  The reason for wanting to create
non-overlapping subregions is to simplify creation of the related masks,
most importantly, the exposure masks.  The \fIoutlimits\fR parameter would
have the values "1 5200", "5201 10400", "10401 15600", and "15601 20800".
The second pair of limits is not specified to obtain the full size along
the second dimension.  Note that the last block will actually be smaller
than 5200 pixels since the survey is less than 20800 pixels.

In each combining step all the images must be specified for the input in
order to compute the full output size but then only those images needed to
produce an output line will be accessed at the same time.  By design this
is roughly 240 images for lines where the strips overlap.  The
non-overlapping blocks can be mosaiced together with this task as a final
step if desired.


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
\fIscale\fR parameter and the zero point shift by the \fIzero\fR
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.

.nf
	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
.fi

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

.nf
	[input|output|overlap] [image section]
.fi

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 value of the
image header keyword specified by the \fIexpname\fR keyword.  As implied
by the parameter name, this is typically the image exposure time since
intensity levels are linear with the exposure time in CCD detectors.
Note that the exposure keyword is also updated in the final image
as the weighted average of the input values.  Thus, if one wants to
use a nonexposure time keyword and keep the exposure time updating
feature the image header keyword syntax is available; i.e. !<keyword>.

Scaling values may be defined as a list of values in a text file.  The file
name is specified by the standard @file syntax.  The list consists of one
value per line.  The order of the list is assumed to be the same as the
order of the input images.  It is a fatal error if the list is incomplete
and a warning if the list appears longer than the number of input images.
Because the scale and zero levels are adjusted only the relative
values are important.

If both an intensity scaling and zero point shift are selected the
zero point is added first and the scaling is done.  This is
important if the scale and offset values are specified by
header keywords or from a file of values.  However,
in the log output the zero values are given as the scale times
the offset hence those numbers would be interpreted as scaling
first and zero offset second.

The image statistics and scale factors are recorded in the log file
unless they are all equal, which is equivalent to no scaling.  The
scale factors are normalized so that the first input image has no
scaling.  This is done because the header of the first input image
is used as the template header for the combined output image.
By scaling to this first image this means that flux related keywords,
such as exposure time and airmass, are representative of the output
(except when the "sum" option is used).

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 \fIsigscale\fR 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 used for the final weighted average, sigma image, and
exposure mask 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:

.nf
	weight_final = weight_initial / (scale * sky)
.fi

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.

When forming the final weighted averages if the sum of the weights of
the non-rejected or excluded pixels is zero then instead of producing
a zero average the unweighted average of the pixels is produced.  Similarly,
in the sigma calculation when the weights of the pixels are all zero
then the sigma is computed as if all pixels have unit weights.

When there are zero weights only the pixels with non-zero weights are
used in computing the output exposure time mask.  Note that the actual
weight values are not used but simply the sum of all exposure times
of pixels from images with non-zero weights is produced.

The purpose of using zero weights is to identify images that are of
poor quality (such as non-photometric or bad seeing) which are then
excluded in the final weighted average or exposure time.  However,
they contribute to the final image when there is no good
quality data but with an output exposure time of zero.

INPUT PIXEL MASKS

A pixel mask is a type of IRAF file having the extension ".pl" or
a FITS extension of "type=mask" which
identifies an integer value with each pixel of the images to which it is
applied.  In future masks may also be stored as special FITS extensions.
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 type of file is compact.
It is 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 \fBmkpattern\fR, edited by \fBimedit\fR, examined by
\fBimexamine\fR, operated upon by \fBimarith\fR, graphed by \fBimplot\fR,
and displayed by \fBdisplay\fR.

To use pixel masks with \fBimcombine\fR 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) or some other
keyword to be specified.  This can be
done with \fBhedit\fR.  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
\fImasktype\fR parameter is set to a value other than "none" or "".  Note that
when it is set to "none", mask information is not used even if it exists for
the image.  The values of \fImasktype\fR which apply masks are "goodvalue",
"badvalue", "novalue", "goodbits", "badbits", and "!<keyword>".  The last choice
allows specifying the keyword whose value is the mask to be used otherwise
the keyword "BPM" is used.

The \fImasktype\fR choices are used in conjunction with the
\fImaskvalue\fR parameter.  When the mask type is "goodvalue" or an
explicit keyword is specified without a mask type, the pixels with mask
values matching the specified value are included in combining and all
others are rejected.  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.

The "novalue" option differs from the others in that there are three
classes of mask values and any output pixel mask will have the three
values 0 for good, 1 for no data, and 2 for bad.  The purpose of this
option is to produce output values from the input values when there are
no good pixels.  This happens when the input images have pixel values
which have been identified as bad (such as saturated) but whose values
can be used, possibly after being replaced or interpolated from nearby
pixels, to produce a value that is either cosmetically reasonable or even
marginally scientifically useful.  Again, this only happens if there
are no good pixels to combine and then the output mask will identify
these pixels with a mask value of 2.  If there is even one good pixel
then only the good data will contribute to the output.  An exposure mask
may be useful in this case when most but not all image pixels have been
eliminated due to things like saturation.

If a mask operation is specified and an image has no mask image associated
with it (the BPM or specified keyword is absent), 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.

When the number of input images exceeds the maximum number of open files
allowed by IRAF and the input images need to be "stacked" then the masks
are also stacked.  The stacking requires all the images to have the same size.


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
\fIlthreshold\fR and \fIhthreshold\fR 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 \fBimedit\fR
or \fBimreplace\fR to set parts of the images to be excluded to some
very low or high magic value.


REJECTION ALGORITHMS

The \fIreject\fR 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.

.in 2
MINMAX
.in 2
A specified fraction of the highest and lowest pixels are rejected.
The fraction is specified as the number of high and low pixels, the
\fInhigh\fR and \fInlow\fR 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,

.nf
	nl = n * nlow/nimages + 0.001 
	nh = n * nhigh/nimages + 0.001 
.fi

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:

.nf
	 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
.fi

.in -2
CCDCLIP
.in 2
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:

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

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 \fBstsdas.wfpc.noisemodel\fR for a way to determine
these values (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 \fIrdnoise\fR 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 \fIgain\fR specifies the gain as either a value or image header
keyword and the parameter \fIsnoise\fR 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 \fIlsigma\fR times this sigma below or \fIhsigma\fR 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 \fIsigscale\fR 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
\fInkeep\fR 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 \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR 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 \fIreject\fR
to select this algorithm, \fImclip\fR to select the median or average for
the center of the clipping, \fInkeep\fR to limit the number of pixels
rejected, the CCD noise parameters \fIrdnoise, gain\fR and \fIsnoise\fR,
\fIlsigma\fR and \fIhsigma\fR to select the clipping thresholds,
and \fIsigscale\fR to set the threshold for making corrections to the sigma
calculation for different image scale factors.

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

.in -2
SIGCLIP
.in 2
The sigma clipping algorithm computes at each output pixel the median or
average excluding the high and low values.  The sigma is then computed
about this estimate (without excluding the low and high values).  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
\fInkeep\fR 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 \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR 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 \fIreject\fR to select
this algorithm, \fImclip\fR to select the median or average for the
center of the clipping, \fInkeep\fR to limit the number of pixels
rejected, \fIlsigma\fR and \fIhsigma\fR to select the
clipping thresholds, and \fIsigscale\fR to set the threshold for
making corrections to the sigma calculation for different image scale
factors.

.in -2
AVSIGCLIP
.in 2
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:

.nf
	sigma(column,line) = sqrt (gain(line) * signal(column,line))
.fi

where the \fIestimated\fR signal is the mean or median (hopefully excluding
any bad pixels) and the gain is the \fIestimated\fR 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
\fInkeep\fR 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 \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR 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 \fIreject\fR to select
this algorithm, \fImclip\fR to select the median or average for the
center of the clipping, \fInkeep\fR to limit the number of pixels
rejected, \fIlsigma\fR and \fIhsigma\fR to select the
clipping thresholds, and \fIsigscale\fR to set the threshold for
making corrections to the sigma calculation for different image scale
factors.

.in -2
PCLIP
.in 2
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 multiplied by the scale factors \fIlsigma\fR and \fIhsigma\fR
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 \fIpclip\fR 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 \fIpclip\fR 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.

After rejection the number of retained pixels is checked against the
\fInkeep\fR 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 \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR parameter only applies to the pixels used
by the clipping rejection algorithm and does not apply to threshold
or bad pixel mask rejection.

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 \fIpclip\fR value of 2 selects the 2nd pixel
above the median (6th) pixel which is the 8th pixel.  A \fIpclip\fR
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 \fIreject\fR to select this
algorithm, \fIpclip\fR to select the percentile pixel, \fInkeep\fR to limit
the number of pixels rejected, and \fIlsigma\fR and \fIhsigma\fR to select
the clipping thresholds.
.in -4

GROW REJECTION

Neighbors of pixels rejected by the rejection algorithms
may also be rejected.  The number of neighbors to be rejected
is specified by the \fIgrow\fR parameter which is a radius in pixels.
If too many pixels are rejected in one of the grown pixels positions
(as defined by the \fInkeep\fR parameter) then the value of that pixel
without growing will be used.

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 as specified by the \fIcombine\fR
parameter.  In all cases if there are no remaining pixels the \fIblank\fR
is produced.  The combining choices are as follows.

.in 2
AVERAGE
.in 2
The weighted average of the remaining pixels is computed.  If no
weighting was specified then a simple, unweighted average is used.
If the sum of the weights of for the accepted pixels is zero then the
unweighted average is output.

.in -2
MEDIAN
.in 2
The median of the remaining pixels is computed.  The median is the
usual mathematical definition where a particular pixel value is produced
for an odd number of pixels and the average of the two central values
is computed for an even number of pixels.

.in -2
SUM
.in 2
The sum of the unrejected pixels is computed.  

.in -2
LMEDIAN
.in 2
The median of the remaining pixels is computed except that for two
pixels the lower value is used.  This is a specialized feature useful
for minimizing the effects of cosmic rays in dithered and/or masked data.

.in -2
QUADRATURE
.in 2
The pixels are combined as

.nf
    sqrt (sum {(wt * sigma)^2}) / sum {wt}
.fi

This is used when the input pixel values represent "sigmas".  This option
is usually a second pass after the input data has been combined.  It is
important that the input is arranged such that the same scaling and
pixel rejections are used.  This means that these cannot be given by
explicit lists and masks and not generated from the data.

.in -2
QUADRATURE
.in 2
The pixels are combined as

.nf
    value = max (0, scaled_pixel_value)
    variance = rdnoise^2 + value / gain + (snoise * value)^2
    output = sqrt (sum {variance * wt^2}) / sum {wt}
.fi

This is used when the variances in the input images can be computed
by the above noise model.  Note that the gain and rdnoise are adjusted
for any scaling applied to the pixel values.

This method has the advantage that the input images are the same as
those used to form a combined image and so all the steps of deriving
scaling and rejecting pixels by some rejection method will be the same.
.in -4

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.
.ih
EXAMPLES
1.  To average and median images without any other features:

.nf
	cl> imcombine obj* avg combine=average reject=none
	cl> imcombine obj* med combine=median reject=none
.fi

2.  To reject cosmic rays:

.nf
	cl> imcombine obs1,obs2 Obs reject=crreject rdnoise=5.1, gain=4.3
.fi

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

.nf
	cl> imcombine @list grid offset="grid 5 65 5 65"
.fi

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

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

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

.nf
	cl> imcombine obj* avsig combine=average reject=avsig \
	>>> scale=exp zero=mode weight=exp  expname=exptime
.fi
.ih
REVISIONS
.ls IMCOMBINE V2.12
A number of enhancements for dealing with large numbers of images were
made.  Also the masktype option "!<keyword>", where <keyword> is a
user specified keyword, was added.

The new parameters "headers", "bpmasks", "rejmasks",  "nrejmasks", and
"expmasks" provide additional types of output.  The old parameters
"rejmask" and "plfile" were removed.  The new "nrejmasks" corresponds
to the old "plfile" and the new "rejmasks" corresponds to the old
"rejmask".

There is a new "combine" type "sum" for summing instead of averaging the
final set of offset, scaled, and weighted pixels.

there is a new parameter "outlimits" to allow output of a subregion of
the full output.  This is useful for raster surveys with large numbers
of images.

Additional keywords may appear in the output headers.

The scaling is now done relative to the first image rather than an
average over the images.  This is done so that flux related keywords
such as exposure time and airmass remain representative.
.le
.ls IMCOMBINE V2.11.2
The grow algorithm was improved to give a 2D growing radius.

An optional output mask file contains the identifications of which pixel
in which input image was rejected or excluded.

The internal calculation type was changed to be the highest precedence
of the input and output types.  Previously it was only the input types.
.le
.ls IMCOMBINE 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.

A combination of ushort and short images now defaults to integer.
.le
.ls IMCOMBINE V2.14
The "masktype" parameter has been generalized to allow both using a
different keyword for the input mask and choosing the mask method.
The "novalue" masktype is new and is useful for maintaining a distinction
between no data and possibly marginally useful or cosmetically useful
data.
.le
.ls IMCOMBINE V2.10.3
The input scalings from an @file or header keyword are now truly
mulitplicative or additive and they are not normalized.  The output
pixel types now include unsigned short integer.
.le
.ls IMCOMBINE 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 \fInkeep\fR parameter was added to allow
controlling the maximum number of pixels to be rejected by the clipping
algorithms.  The \fIsnoise\fR parameter was added to include a sensitivity
or scale noise component to the noise model.  Errors will now delete
the output images.
.le
.ls IMCOMBINE V2.10
This task was greatly revised to provide many new features.  These features
are:

.nf
    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
.fi
.le
.ih
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.
.ih
SEE ALSO
ccdred.combine mscred.combine onedspec.scombine, wpfc.noisemodel,
obsolete.ocombine
.endhelp