FUNCTION xrt_despike, image_in, sens=sens, spike_map=ss, hist=hist, $
                      verbose=verbose, quiet=quiet, qabort=qabort,  $
                      qstop=qstop

; =========================================================================
;+
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;       XRT_DESPIKE
;
; CATEGORY:
;       Data calibration
;
; PURPOSE:
;       To smooth out "bad" pixels due to cosmic ray hits. Generally, 
;       this is intended as a cosmetic fix, not as a scientific
;       recalibration of the spiked pixels. However, unspiked pixels
;       are unaffected, and the SPIKE_MAP can be used to identify the
;       corrected pixels. See "Notes" for more information.
;
; CALLING SEQUENCE:
;       Result = XRT_DESPIKE(image_in [,sens=sens] [,spike_map=spike_map]
;                            [,hist=hist] [,/verbose] [,/quiet] [,/qstop] 
;                            [,qabort=qabort] )
;
; INPUTS:
;
;       IMAGE_IN  - [Mandatory] (2-dim number array, [Nx,Ny])
;                   The input image array.
;
; KEYWORDS:
;
;       SENS      - [Optional] (number scalar)
;                   Sensitivity to spots. If a pixel has a value greater
;                   than (SENS * Mean(neighbors)), then it is identified
;                   as a spike. Values of 1.1 to 1.30 are reasonable.
;                   Default value = 1.1.
;       /VERBOSE  - [Optional] (Boolean) If set, print out extra
;                   information. Overrides "/quiet" (see Note #3).
;       /QUIET    - [Optional] (Boolean) If set, suppress messages
;                   (see Note #3).
;       /QSTOP    - [Optional] (Boolean) For debugging.
;
; OUTPUTS:
;
;       Return    - [Mandatory] (2-dim number array, [Nx,Ny])
;                   The output image array, with the spikes reduced.
;                   A spike value is replaced with the median of the 16
;                   pixels marking a 5x5 square around the spike.
;       SPIKE_MAP - [Optional] (Long array, [Nspikes])
;                   The addresses of the despiked pixels. These are the
;                   only pixels altered by this routine.
;       HIST      - [Optional] (String scalar)
;                   Brief report on how many spiked pixels were found.
;                   Called "hist" because it fits the format for an
;                   index.HISTORY entry. (This keyword is used by 
;                   <xrt_prep.pro>, for example.)
;       QABORT    - [Optional] (Boolean) Indicates that the program
;                   exited gracefully without completing. (Might be
;                   useful for calling programs.)
;                   0: Program ran to completion.
;                   1: Program aborted before completion.
;
; EXAMPLES:
;
;       Basic usage:
;       IDL> image_out = xrt_despike(image_in)
;
;       Find out which pixels were altered:
;       IDL> image_out = xrt_despike(image_in, spike_map=spike_map)
;       IDL> help, image_in[spike_map]
;
;
; COMMON BLOCKS:
;       None.
;
; NOTES:
;
;       1) Spiked pixels are identified if they are too much larger
;          than their average neighbor. (See SENS.) Only such pixels are
;          altered by this routine. They are replaced with the median
;          of a slightly-removed ring of pixels. (See "Return".)
;
;       2) Multiple passes may be required to affect heavy image
;          spiking. Generally, one pass does fine for lightly spiked
;          images, at SENS=1.1. A second pass should handle most of the
;          heavily spiked images. Picky people may need to experiment a 
;          bit with the number of passes and the SENS level.
;
;       3) There are three levels of verbosity.
;          a) "verbose" = highest priority. All errors and messages are
;                         displayed. ("if q_vb")
;          b) "quiet" = lower priority. No errors or messages are
;                       displayed. ("if q_qt")
;          c) neither = lowest priority. All errors and some messages are
;                       displayed. ("if not q_qt")
;
; MODIFICATION HISTORY:
;
progver = 'v2007.May.18' ; --- (P.Jibben, M.Weber) Written. Jibben derived
;                              the cleaning algorithms based upon the 
;                              ANA <browser> software and "The Handbook of
;                              Astronomical Image Processing" by Berry &
;                              Burnell. Weber implemented a faster version
;                              of the code.
;
;-
; =========================================================================

;=== Initial setup =======================================

  ;=== Set Booleans which control print statements.
  ;=== Keyword "verbose" overrules "quiet".
  q_vb = keyword_set(verbose)
  q_qt = keyword_set(quiet) and (q_vb eq 0)

  ;=== Initialize program constants.
  prognam='XRT_DESPIKE'
  prognul='           '
  def_sens = 1.1

  ;=== Announce version of the program.
  if q_vb then box_message, prognam+': Running ' + progver + '.'
  if q_vb then print, prognam+': Performing initial setup...'

  ;=== Set some keyword Booleans.
  qstop   = keyword_set(qstop)
  qabort  = 0B

  if q_vb then print, prognam+'                          ...OK'

;=== Check inputs ========================================

  if q_vb then print, prognam+': Checking inputs...'

  ;=== Check for number of inputs.
  if (n_params() ne 1) then begin
    if (not q_qt) then box_message, prognam+': Incorrect number of ' $
                                    + 'parameters. Aborting.'
    qabort = 1B
    return, 0
  endif

  ;=== Check for image dimensions.
  imsiz = size(image_in)
  if (imsiz[0] ne 2) then begin
    if (not q_qt) then box_message, prognam+': Input IMAGE_IN must ' $
                                        + 'be a 2D image. Aborting.'
    qabort = 1B
    return, 0
  endif

  ;=== Check SENS value.
  if keyword_set(sens) then begin

    if (n_elements(sens) ne 1) then begin
      if (not q_qt) then box_message, [prognam+': Input SENS must ' $
        + 'be a scalar.', prognul+'  Using the default value = ' $
        + strcompress(def_sens,/rem) + '.']
      sens = def_sens
    endif

    if (not is_number(sens)) then begin
      if (not q_qt) then box_message, [prognam+': Input SENS must ' $
        + 'be a number.', prognul+'  Using the default value = ' $
        + strcompress(def_sens,/rem) + '.']
      sens = def_sens
    endif

    if (sens lt 1.0) then begin
      if (not q_qt) then box_message, [prognam+': Input SENS must ' $
        + 'be a number > 1.0.', prognul+'  Using the default value = ' $
        + strcompress(def_sens,/rem) + '.']
      sens = def_sens
    endif

  endif else begin
    sens = def_sens
  endelse

  if q_vb then print, prognul+'                 ...OK'

;=== Find and remove spikes ==============================

  if q_vb then print, prognam+': Fixing spikes...'

  nx = n_elements(image_in[*,0])
  ny = n_elements(image_in[0,*])

  ;=== Find avg for nearest neighbors.
  kernel = [[1,1,1],[1,0,1],[1,1,1]] / 8.
  img_avg = convol(image_in, kernel)

  ;=== Compare pix/avg versus threshold,
  ;=== for pixels not near the edge. 
  ss = where( (image_in[2:nx-3,2:ny-3] / img_avg[2:nx-3,2:ny-3]) $
             ge sens, count)

  if q_vb then print, prognam+': Found '+strcompress(count,/rem) $
                      +' pixels to be fixed.'

  ;=== Found some spikes? Then proceed.
  if (count ge 1) then begin

    ;=== Have to convert ss to be relative to entire image.
    aa = bytarr(nx-4,ny-4)
    aa[ss] = 1
    bb = bytarr(nx,ny)
    bb[2:nx-3,2:ny-3] = aa
    ss = where(bb eq 1, count)

    ;=== Prepare an "address" filter for +2 pixels around spikes.
    filter = [-2*nx-2, -2*nx-1, -2*nx, -2*nx+1, -2*nx+2, -nx-2, -nx+2, $
              -2, 2, nx-2, nx+2, 2*nx-2, 2*nx-1, 2*nx, 2*nx+1, 2*nx+2   ]

    ;=== spike_ring is the filter applied to each spike address.
    spike_ring = lonarr(16,count)
    for ii = 0,15 do spike_ring[ii,*] = ss + filter[ii]

    ;=== For each spike, find median of its ring.
    ;=== This bit lets us do them all simultaneously.
    spike_med = median(image_in[spike_ring], dim=1, /even)

    ;=== Replace spikes with medians.
    image_out     = image_in
    image_out[ss] = spike_med

  ;=== Else if didn't find any spikes.
  endif else image_out = image_in

  if q_vb then print, prognul+'               ...OK'

;=== Finish ==============================================

  hist = prognam+': Despiked '+strcompress(count,/rem)+' pixels, at a ' $
         + 'sensitivity of '+strcompress(sens,/rem)+'.'

  return, image_out

END ;======================================================================