CR_REJECT
Name
CR_REJECT
Purpose
General, iterative cosmic ray rejection using two or more input images.
Explanation
Uses a noise model input by the user, rather than
determining noise empirically from the images themselves.
The image returned has the combined exposure time of all the input
images, unless the bias flag is set, in which case the mean is
returned. This image is computed by summation (or taking mean)
regardless of loop and initialization options (see below).
Calling Sequence
cr_reject, input_cube, rd_noise_dn, dark_dn, gain, mult_noise, $
combined_image, combined_npix, combined_noise
MODIFIED ARGUMENT:
input_cube - Cube in which each plane is an input image.
If the noise model is used (rd_noise_dn, dark_dn,
gain), then input_cube must be in units of DN.
If an input noise cube is supplied (rd_noise_dn
<0), then the units of input_cube and noise_cube
merely need to be consistent.
This array is used as a buffer and its contents
are not guaranteed on output (although for now,
weighting=0 with /restore_sky should give you back
your input unaltered).
Input Arguments
rd_noise_dn - Read noise per pixel. Units are DN.
If negative, then the user supplies an error cube
via the keyword noise_cube. In the latter case,
mult_noise still applies, since it is basically a fudge.
dark_dn - Dark rate in DN per pixel per s. This can be a scalar,
or it can be a dark image divided by the exposure
time.
gain - Electrons per DN.
mult_noise - Coefficient for multiplicative noise term -- helps
account for differing PSFs or subpixel image shifts.
Input Keywords
exptime - If the images have different exposure times, pass
them in a vector. You can leave this off for
frames with the same exposure time, but dark counts
won't be treated correctly.
verbose - If set, lots of output.
nsig - Rejection limit in units of pixel-to-pixel noise
(sigma) on each input image. Can be specified as
an array, in which case the dimension gives the
maximum number of iterations to run. (Default =
[8, 6, 4]
dilation - With dfactor, provides functionality similar to the
expansion of the CR with pfactor and radius in STSDAS
crrej. Dilate gives the size of the border to be
tested around each initially detected CR pixel.
E.g., dilate=1 searches a 9 X 9 area centered on the
original pixel. If dfactor is set, the default is 1.
dfactor - See dilation. This parameter is equivalent to pfactor
in STSDAS crrej. The current threshold for rejection
is multiplied by this factor when doing the search
with the dilated mask. If dilation is set, the default
for this parameter is 0.5.
bias - Set if combining biases (divides through by number
of images at end, since exposure time is 0).
tracking_set - Subscripts of pixels to be followed through the
computation.
noskyadjust - Sky not to be subtracted before rejection tests. Default
is to do the subtraction.
xmedsky - Flag. If set, the sky is computed as a 1-d array
which is a column-by-column median. This is intended
for STIS slitless spectra. If sky adjustment is
disabled, this keyword has no effect.
input_mask - Mask cube input by the user. Should be byte data
because it's boolean. 1 means use the pixel,
and 0 means reject the pixel - these rejections
are in addition to those done by the CR rejection
algorithm as such.
The following keywords control how the current guess at a CR-free
"check image" is recomputed on each iteration:
median_loop - If set, the check image for each iteration is
the pixel-by-pixel median. THE MEAN IS
RETURNED in combined_image even if you set
this option. (Default is mean_loop.)
minimum_loop - If set, the check image for each iteration is
the pixel-by-pixel minimum. THE MEAN IS
RETURNED in combined_image even if you set
this option. (Default is mean_loop.)
mean_loop - If set, the check image for each iteration is
the pixel-by-pixel mean. (Same as the default.)
noclearmask - By default, the mask of CR flags is reset before
every iteration, and a pixel that has been
rejected has a chance to get back in the game
if the average migrates toward its value. If this
keyword is set, then any rejected pixel stays
rejected in subsequent iterations. Note that what
stsdas.hst_calib.wfpc.crrej does is the same
as the default. For this procedure, the default
was NOT to clear the flags, until 20 Oct. 1997.
restore_sky - Flag. If set, the routine adds the sky back into
input_cube before returning. Works only if
weighting=0.
null_value - Value to be used for output pixels to which no
input pixels contribute. Default=0
weighting - Selects weighting scheme in final image
combination:
0 (default) - Poissonian weighting - co-add
detected DN from non-CR pixels. (Pixel-by-
pixel scaling up to total exposure time,
for pixels in stack where some rejected.)
Equivalent to exptime weighting of rates.
1 or more - Sky and read noise weighting of rates.
Computed as weighted average of DN rates,
with total exp time multiplied back in
afterward.
In all cases, the image is returned as a sum in
DN with the total exposure time of the image
stack, and with total sky added back in.
The following keywords allow the initial guess at a CR-free "check
image" to be of a different kind from the iterative guesses:
init_med - If set, the initial check image is
the pixel-by-pixel median. (Not permitted if
input_cube has fewer than 3 planes; default is minimum.)
init_mean - If set, the initial check image is
the pixel-by-pixel mean. (Default is minimum.)
init_min - If set, the initial check image is
the pixel-by-pixel minimum. (Same as the default.)
Output Arguments
combined_image - Mean image with CRs removed.
combined_npix - Byte (or integer) image of same dimensions as
combined_image, with each element containing
the number of non-CR stacked pixels that
went into the result.
combined_noise - Noise in combined image according to noise model
or supplied noise cube.
Output Keywords
mask_cube - CR masks for each input image. 1 means
good pixel; 0 means CR pixel.
skyvals - Sky value array. For an image cube with N planes,
this array is fltarr(N) if the sky is a scalar per
image plane, and fltarr(XDIM, N), is the XMEDSKY
is set.
INPUT/OUTPUT KEYWORD:
noise_cube - Estimated noise in each pixel of input_cube as
returned (if rd_noise_dn ge 0), or input noise
per pixel of image_cube (if rd_noise_dn lt 0).
skybox - X0, X1, Y0, Y1 bounds of image section used
to compute sky. If supplied by user, this
region is used. If not supplied, the
image bounds are returned. This parameter might
be used (for instance) if the imaging area
doesn't include the whole chip.
COMMON BLOCKS: none
SIDE EFFECTS: none
Method
COMPARISON WITH STSDAS
Cr_reject emulates the crrej routine in stsdas.hst_calib.wfpc.
The two routines have been verified to give identical results
(except for some pixels along the edge of the image) under the
following conditions:
no sky adjustment
no dilation of CRs
initialization of trial image with minimum
taking mean for each trial image after first (no choice
in crrej)
Dilation introduces a difference between crrej and this routine
around the very edge of the image, because the IDL mask
manipulation routines don't handle the edge the same way as crrej
does. Away from the edge, crrej and cr_reject are the same with
respect to dilation.
The main difference between crrej and cr_reject is in the sky
computation. Cr_reject does a DAOPHOT I sky computation on a
subset of pixels grabbed from the image, whereas crrej searches
for a histogram mode.
REMARKS ON USAGE
The default is that the initial guess at a CR-free image is the
pixel-by-pixel minimum of all the input images. The pixels
cut from each component image are the ones more than nsig(0) sigma
from this minimum image. The next iteration is based on the
mean of the cleaned-up component images, and the cut is taken
at nsig(1) sigma. The next iteration is also based on the mean with
the cut taken at nsig(2) sigma.
The user can specify an arbitrary sequence of sigma cuts, e.g.,
nsig=[6,2] or nsig=[10,9,8,7]. The user can also specify that
the initial guess is the median (/init_med) rather than the
minimum, or even the mean. The iterated cleaned_up images after
the first guess can be computed as the mean or the median
(/mean_loop or /median_loop). The minimum_loop option is also
specified, but this is a trivial case, and you wouldn't want
to use it except perhaps for testing.
The routine takes into account exposure time if you want it to,
i.e., if the pieces of the CR-split aren't exactly the same.
For bias frames (exposure time 0), set /bias to return the mean
rather than the total of the cleaned-up component images.
The crrej pfactor and radius to propagate the detected CRs
outward from their initial locations have been implemented
in slightly different form using the IDL DILATE function.
It is possible to end up with output pixels to which no valid
input pixels contribute. These end up with the value
NULL_VALUE, and the corresponding pixels of combined_npix are
also returned as 0. This result can occur when the pixel is
very noisy across the whole image stack, i.e., if all the
values are, at any step of the process, far from the stack
average at that position even after rejecting the real
outliers. Because pixels are flagged symmetrically N sigma
above and below the current combined image (see code), all
the pixels at a given position can end up getting flagged.
(At least, that's how I think it happens.)
Modification History
5 Mar. 1997 - Written. R. S. Hill
14 Mar. 1997 - Changed to masking approach to keep better
statistics and return CR-affected pixels to user.
Option to track subset of pixels added.
Dilation of initially detected CRs added.
Other small changes. RSH
17 Mar. 1997 - Arglist and treatment of exposure times fiddled
to mesh better with stis_cr. RSH
25 Mar. 1997 - Fixed bug if dilation finds nothing. RSH
4 Apr. 1997 - Changed name to cr_reject. RSH
15 Apr. 1997 - Restyled with emacs, nothing else done. RSH
18 Jun. 1997 - Input noise cube allowed. RSH
19 Jun. 1997 - Multiplicative noise deleted from final error. RSH
20 Jun. 1997 - Fixed error in using input noise cube. RSH
12 July 1997 - Sky adjustment option. RSH
27 Aug. 1997 - Dilation kernel made round, not square, and
floating-point radius allowed. RSH
16 Sep. 1997 - Clearmask added. Intermediate as well as final
mean is exptime weighted. RSH
17 Sep. 1997 - Redundant zeroes around dilation kernel trimmed. RSH
1 Oct. 1997 - Bugfix in preceding. RSH
21 Oct. 1997 - Clearmask changed to noclearmask. Bug in robust
array division fixed (misplaced parens). Sky as
a function of X (option). RSH
30 Jan. 1998 - Restore_sky keyword added. RSH
5 Feb. 1998 - Quick help corrected and updated. RSH
6 Feb. 1998 - Fixed bug in execution sequence for tracking_set
option. RSH
18 Mar. 1998 - Eliminated confusing maxiter spec. Added
null_value keyword. RSH
15 May 1998 - Input_mask keyword. RSH
27 May 1998 - Initialization of minimum image corrected. NRC, RSH
9 June 1998 - Input mask cube processing corrected. RSH
21 Sep. 1998 - Weighting keyword added. RSH
7 Oct. 1998 - Fixed bug in input_mask processing (introduced
in preceding update). Input_mask passed to
skyadj_cube. RSH
5 Mar. 1999 - Force init_min for 2 planes. RSH
1 Oct. 1999 - Make sure weighting=1 not given with noise cube. RSH
1 Dec. 1999 - Corrections to doc; restore_sky needs weighting=0. RSH
17 Mar. 2000 - SKYBOX added. RSH