>  Docs Center  >  Libraries  >  ASTROLIB  >  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).

Calling Sequence

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

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
    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
    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
    null_value - Value to be used for output pixels to which no
                    input pixels contribute. Default=0
    weighting - Selects weighting scheme in final image
                    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
                    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.
    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.


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

© NV5 Geospatial Solutions, Inc. |  Legal
   Contact Us