cr_reject.pro

basic median filter simulation

PRO cr_reject, input_cube, rd_noise_dn, dark_dn, gain, mult_noise, $
               combined_image, combined_noise, combined_npix, $
               MASK_CUBE=mask_cube, NOISE_CUBE=noise_cube, $
               NSIG=nsig, MEDIAN_LOOP=median_loop, MEAN_LOOP=mean_loop, $
               MINIMUM_LOOP=minimum_loop, INIT_MED=init_med, $
               INIT_MIN=init_min, INIT_MEAN=init_mean, EXPTIME=exptime,$
               BIAS=bias, VERBOSE=verbose, $
               TRACKING_SET=tracking_set, DILATION=dilation, DFACTOR=dfactor, $
               NOSKYADJUST=noskyadjust,NOCLEARMASK=noclearmask, $
               XMEDSKY=xmedsky, RESTORE_SKY=restore_sky, $
               SKYVALS=skyvals, NULL_VALUE=null_value, $
               INPUT_MASK=input_mask, WEIGHTING=weighting, SKYBOX=skybox
;+
; 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
;-
on_error,0
IF n_params(0) LT 6 THEN BEGIN
    print,'CALLING SEQUENCE:  cr_reject, input_cube, rd_noise_dn, $'
    print,'   dark_dn, gain, mult_noise, combined_image, combined_noise, $'
    print,'   combined_npix'
    print,'KEYWORD PARAMETERS:  nsig, exptime, bias, verbose,'