Automatic Cosmic Ray Extraction
  This program will attempt to identify and remove Cosmic Ray strikes from
    an image.  This program was developed and tested on HST PC data prior
    to the refurbishment mission.  It may work for other types of data but
    it is as of now untested elsewhere.
  The simplest usage is the single pass mode where the same parameters are
    used for the entire image.  First, the image is smoothed with a median
    filter using a box filter size given by width.  This smoothed image is
    then subtracted from the original image.  A robust mean of a portion of
    the image is calculated and subtracted from the image though this mean
    should be near zero.  Any pixel see to deviate by THRESH standard
    deviations from this average is marked for removal and replaced by its
    corresponding value in the smoothed image.
  This initial step works very well on the sky.  I've found that THRESH=3
    WIDTH=7 work pretty well on all but the largest CRS's.  Using a value
    for width less than 7 seems to leave residual "rings" of hot pixels from
    around the edges of a strike.
  The draw back to these parameters is that it is much too agressive in and
    around actual objects in the frame.  The cores of the PSF will be removed
    and numerous pixels will be tagged in the wings of the PSF.
  To get around this problem, use the EXCLUDE keyword.  This is a 5xN array
    containing circular regions to scan with different parameters.  The array
    holds N such regions.  For each region you must specify the following:
       (0,n) - x location of region
       (1,n) - y location of region
       (2,n) - radius of region
       (3,n) - sigma threshold to use in this region
       (4,n) - width of median smoothing, (no smoothing if set to 0).
    The procedure used is to first restore all the pixels in the region to
    their original values (in case they were changed in the first step).
    If the width is set to zero, nothing more is done.  If the width is a
    meaningful value, then the original image is smoothed with that width
    and the region is scanned for deviant pixels again and replaces any
  The effects of each region are cumulative on the image and done in order
    they appear in the array.  So any final steps of restoring small locations
    should be done last.  Also, the step of smoothing the image is very cpu
    intensive.  It will run much faster if you can group all regions with
    similar smoothing values together.
  In practice, a smaller width (~3) and higher thresh (~4) seems to work in
    the wings of the PSF, but it will still take out the core.  So, you
    need to specify two zones, one for cleaning that is nearly the size
    of the outer fringes of the PSF and one for pretecting the image that
    is smaller and centered on the core.  If you happen to get a strike
    near the core of the PSF, this routine won't help and you're probably
    screwed anyway.
  CCD data processing
  pro acre,dirty_im,clean_im,thresh,width, $
              BLFINAL=blfinal,BLMASK=blmask, $
  dirty_im - Original input image to be cleaned.
  thresh   - Deviation threshold, in sigma, from background to cause
                pixel to be fixed.
  width    - Median smoothing width to get local background reference.
  BLFINAL  - Flag, if true, brings up ITOOL to blink between the original
                and cleaned images.
  BLMASK   - Flag, if true, brings up ITOOL to blink between the original
                and the mask showing pixels that are being replaced.
  EXCLUDE  - Array that controls special extraction behavior in select
                regions of the image.  See DESCRIPTION for details.
  VERBOSE  - Flag, if true, generates a wordy output of progress and
                action as routine progresses.
  clean_im - Final cleaned up image.
  MASK     - Return of the mask image.
  94/04/05 - Written by Marc W. Buie, Lowell Observatory