CR_REJECT

     General, iterative cosmic ray rejection using two or more input images.

     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).

     cr_reject, input_cube, rd_noise_dn, dark_dn, gain, mult_noise, $
        combined_image, combined_npix, combined_noise

     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).

     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.

     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.

     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.

     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

     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.)

      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

[Home Page] [Software, Documentation] [IDL Documentation] [Quick Reference] [Feedback]