; Copyright (c) 1998-2010 C.J. Jacobsen  All rights reserved
;+
;NAME:
;   DESPIKE.PRO
;
;LAST CHANGED: -----------------------------------  18-Jul-10
;
; PURPOSE:
; This function replace these higher-than-threshold pixels with the
; average of their neighbors
;
; CATEGORY:
; Image processing - called in zstack_analyze
;
; CALLING SEQUENCE:
; Result = DESPIKE(image, [median_percent=median_percent,$
;                leefilt_percent=leefilt_percent,help=help])
;
; INPUTS:
; Image   2-d array
; 
; KEYWORDS:
; median_percent   sets median_percent threshold
; leefilt_percent sets leefilt_percent threshold (default =  50%)
; help            prints help file
;
; OUTPUTS:
; Returns despiked image
; 
; COMMON BLOCKS: none
;
; MODIFICATION HISTORY:
; (28-Nov-03 cjj) last modified (file extracted from IDL-X1.ZIP of 26-Feb-2009)
; (20-Jul-10 aph) AXIS standard header added
;-

FUNCTION despike,image,$
                 median_percent=median_percent,$
                 leefilt_percent=leefilt_percent,$
                 help=help
  
  IF (keyword_set(help) OR (n_params() EQ 0)) THEN BEGIN
      print,'despiked_image=despike(image,median_percent=,leefilt_percent=,threshold=)'
      print,'  Finds a threshold value by one of two methods:'
      print,'    1. Calculate the LEEFILT version of the image, and find'
      print,'       pixels that are leefilt_percent above the maximum'
      print,'       (default, with leefilt_percent=50)'
      print,'    2. Calculate the MEDIAN(,3) version of the image, and find'
      print,'       pixels that are median_percent above the maximum'
      print,'  Then replace these higher-than-threshold pixels with the'
      print,'  average of their neighbors.'
      return,0.
  ENDIF
  
  svec = size(image)
  IF (svec(0) NE 2) THEN BEGIN
      message,'despike(image) wants a 2D array'
      return, 0.
  ENDIF
  nx=svec(1)
  ny=svec(2)
  
  IF keyword_set(median_percent) THEN BEGIN
      median_max = max(median(image,3))
      threshold = (1.+0.01*median_percent)*median_max
  ENDIF ELSE BEGIN
      svec = size(leefilt_percent)
      IF (svec(svec(0)+1) EQ 0) THEN BEGIN
          leefilt_percent = 50.
      ENDIF
      leefilt_max = max(leefilt(image))
      threshold = (1.+0.01*leefilt_percent)*leefilt_max
  ENDELSE
  spikes = where(image GT threshold, n_spikes)
  resulting_image = image
  IF (n_spikes NE 0) THEN BEGIN
      FOR i_spike=0,(n_spikes-1) DO BEGIN
          ix = spikes(i_spike) mod nx
          iy = spikes(i_spike) / nx
          IF (ix EQ 0) THEN BEGIN
              ix1 = 1
              ix2 = 2
          ENDIF ELSE IF (ix EQ (nx-1)) THEN BEGIN
              ix1 = nx-2
              ix2 = nx-3
          ENDIF ELSE BEGIN
              ix1 = ix-1
              ix2 = ix+1
          ENDELSE
          
          IF (iy EQ 0) THEN BEGIN
              iy1 = 1
              iy2 = 2
          ENDIF ELSE IF (iy EQ (ny-1)) THEN BEGIN
              iy1 = ny-2
              iy2 = ny-3
          ENDIF ELSE BEGIN
              iy1 = iy+1
              iy2 = iy-1
          ENDELSE
          
          resulting_image(ix,iy) = (0.25*(image(ix1,iy)+image(ix2,iy)+$
                                          image(ix,iy1)+image(ix,iy2)))
      ENDFOR
  ENDIF
  
  return, resulting_image
END


