;+
; $Id: despike_gen.pro,v 1.2 2007/02/07 20:42:40 thompson Exp $
;
; NAME:
;        despike_gen
; PURPOSE:
;        Remove (cosmic ray) spikes from an image.
;
;        The basic idea is to median filter the image.  If the difference
;        between the image and its median exceeds a threshold, then that
;        image pixel is replaced by the median.
;
;        This program features image contents and image noise dependent
;        thresholds, plus a more generalized concept of the median which
;        proves to be particularly powerful if cosmic ray hits cause pixel
;        bleeding.
;        The program can be run repeatedly on its output, which allows
;        removal of artifacts significantly larger than one pixel.
;
;        The optimal choice of parameters depends on the instrument, e.g.,
;        on the telescope's point-spread function, the instrument gain
;        (DN per photon), and others.  Specific instrument flags set
;        useful defaults derived from previous experience.
;
;        IMPORTANT: despike_gen should be run on "near-raw" data.  It can
;                   be dark subtracted and flat-fielded, but it should not
;                   have gone through massaging such as (non-linear)
;                   rescaling, smoothing, edge-enhancement, etc.
;                   In addition, the image SHOULD NOT BE COMPRESSED.
;                   The algorithm depends on a linear intensity scale to
;                   predict the photon noise, and on the image more or less
;                   reflecting the instrument point-spread function.
; CATEGORY:
; CALLING SEQUENCE:
;        out = despike_gen(im)
; INPUTS:
;        im = input image
; KEYWORDS (INPUT)
;        Note: the threshold is the maximum of 3 components, controlled by the
;              3 keyword parameters (below):
;              threshold = t_const > t_noise * sqrt(im) > t_fluct * med_fluct
;                          with: med_fluct = median(abs(im-median(im)))
;              The optimal parameter values depend on the instrument, and
;              (ideally) not on the particular image.
;        tconst = scalar.  Constant threshold.
;                  Default: t_const = 3
;        tnoise = scalar.  Scale factor for photon noise dependent threshold.
;                  Default: t_noise = 3
;                  One may want to set t_noise to at least 3 * #DN/photon, or
;                  higher, especially if /low3 is used.
;        tfluct = scalar.  Scale factor for signal fluctuation dependent thr.
;                  Default: t_fluct = 0  (no threshold)
;                  A non-zero value of t_fluct can help reduce the number of
;                  real features being changed, in particular in heavily
;                  undersampled data where t_noise has to be set rather low.
;        /low3 = flag.  If set, the 3x3 median filter is replaced by a filter
;                that picks the 3rd lowest value of the 3x3 pixel neighborhood
;                (median picks the 5th lowest).  This is particularly useful
;                for removing extended artefacts (combined with running
;                despike_gen repeatedly), or in the presence of large numbers
;                of cosmic ray hits.
;                /low3 causes the program to run much slower however.
;        /eit = flag.  Default settings that seem to work with EIT data.
;               These initial settings were only tested with half resol. data
;        /trace = flag.  Default settings that seem to work with TRACE data.
; OUTPUTS:
;        out = despiked image
; KEYWORDS (OUTPUT):
;        events = vector of addresses of the spike pixels
;        threshold = threshold used (image)
; COMMON BLOCKS:
;        None.
; SIDE EFFECTS:
;        Calls "med3x3gen" function
; RESTRICTIONS:
; PROCEDURE:
;        Spikes are those pixels whose intensity exceeds the 3x3 median
;        by a threshold.  The threshold can be
;        - a constant tconst, or
;        - tnoise times the square root of the intensity, or
;        - tfluct times the median fluctuation in a 3x3 neighborhood, or
;        the largest of any combination of the above.
;        A modified median can be used which picks the 3rd lowest pixel
;        in a 3x3 neighborhood (/low3).
; MODIFICATION HISTORY:
;        JPW, 28-apr-98
;
; $Log: despike_gen.pro,v $
; Revision 1.2  2007/02/07 20:42:40  thompson
; WTT - Fixed bug of WHERE function not finding any matches
;
; Revision 1.1  2007/01/25 23:23:53  euvi
; initial 1998 version
;
;-

function despike_gen,im,tconst=t_c,tnoise=t_n,tfluct=t_f,low3=low3, $
   eit=eit,trace=trace,events=eve,threshold=thr

; set proper defaults
if keyword_set(eit) then begin
   if n_elements(t_c) eq 0 then t_c = 3
   if n_elements(t_f) eq 0 then t_f = 3
   if n_elements(t_n) eq 0 then $
      if keyword_set(low3) then t_n = 6 else t_n = 3
endif
if keyword_set(trace) then begin
   if n_elements(t_c) eq 0 then t_c = 3
   if n_elements(t_f) eq 0 then t_f = 0
   if n_elements(t_n) eq 0 then t_n = 3
   if n_elements(low3) eq 0 then low3 = 1
endif
if n_elements(t_c) eq 0 then t_c = 3
if n_elements(t_f) eq 0 then t_f = 0
if n_elements(t_n) eq 0 then t_n = 3

; If /low3 set, use the 3rd lowest = 7th highest in the 3x3 neighborhood
; instead of the median (5th highest)
; ame = amedian(im,3)
if keyword_set(low3) then pick=7 else pick=5
ame = med3x3gen(im,pick=pick)

; calculate median fluctuation
; flu = amedian(abs(im-ame),3)
flu = med3x3gen(abs(im-ame),pick=5)

; calculate (photon) noise
noi = sqrt(im > 0)

thr = (t_c > t_f*flu) > t_n*noi

eve = where(im-ame gt thr, count)

out = im
if count gt 0 then out(eve) = ame(eve)

return,out

end