;+  
;
; Laboratory for Atmospheric and Space Physics
; University of Colorado, Boulder, Colorado, USA
;
; FILENAME:
;   remove_spikes.pro
;
; AUTHOR:
;   Lon Riesberg
;
; DATE: November 15, 2005 
;
; PURPOSE:
;   Removes spikes from calibrated and corrected level 1A images.  
;
; INPUT PARAMETERS:
;   images - array of image structures that need to be cleaned of spikes.
;
; OUTPUT PARAMETERS:
;   NONE
;
; RETURN VALUE:
;   images - array of image structures where spikes have been removed.
;  
; BACKGROUND:
;
; ALGORITHM:
;      1. A sliding window checks a 'cluster' of pixels for potential spikes.
;         A potential spike is identified as any pixels that's at least a
;         default number of sigmas difference from the cluster mean.
;      2. Pixels that have been flagged as potential spikes must be flagged
;         as potential spikes in every cluster it participates in before
;         flagged as a true spike. 
;      3. The values of any pixels flagged as spikes get set to NAN.
;
; REFERENCES:
;   CIPS Algorithm Theoretical Basis Document
;   Conversations with CIPS Science and Engineering Teams(Bill McClintock) 
; 
; NOTES:
; 
; CONSTRAINTS:
; 
; OTHER SERVERS/MODULES USED:
;
;  RELATED MODULES / CLASSES:
;     NONE
;  
;  SUPPORTING DATABASE TABLES OR FILES:
;
;  USAGE EXAMPLE:
;	spike_corrected_images = remove_spikes(images)      
; 

function remove_spikes, images

   ;
   ; INIT
   ;
   ; size of pixel clusters in each direction from any given pixel
   ; (when cluster_size_x is 2 and cluster_size_y is 2, the total
   ; cluster size is 5 x 5 pixels)
   cluster_size_x = 2
   cluster_size_y = 2

   ; number of sigmas difference from cluster mean to identify a 
   ; given pixel as being a spike relative to other pixels in a cluster
   spike_defn = 5 

   ;
   ; FOR EACH IMAGE, IDENTIFY SPIKES
   ;
   for i=0, n_elements(images)-1 do begin

      ; get dimensions of current image
      dim = size(*(images[i].image), /dimensions)
      xdim = dim[0]
      ydim = dim[1]
   
      ; array to flag potential spikes
      ; each pixel participates in a number of clusters... num_times_checked
      ; is an array that keeps track of the number of times each pixel
      ; is checked as being a spike which essentially counts the number of 
      ; clusters each pixel participates in
      num_times_checked = make_array(dim, /byte, value = 0)
      num_times_flagged_as_spike = make_array(dim, /byte, value = 0)
   
      ; for each pixel, check if current pixel is a spike relative to its neighbors
      for x=cluster_size_x, xdim-cluster_size_x-1 do begin
         for y=cluster_size_y, ydim-cluster_size_y-1 do begin

            ; extract the cluster of pixels that x, y belongs with
            start_x = x - cluster_size_x         
            end_x = x + cluster_size_x     
            start_y = y - cluster_size_y
            end_y = y + cluster_size_y 
            pixel_cluster = (*(images[i].image))[start_x:end_x, start_y:end_y]
   
            ; update num_times_checked array
            num_times_checked[start_x:end_x, start_y:end_y] += 1

            ; determine if any pixels in the cluster are spikes
            wild_indices = where(not_wild(pixel_cluster, spike_defn) eq 0, num_wild)

            if num_wild ne 0 then begin                 
               for j=0, n_elements(wild_indices)-1 do begin
                  wild_x = (array_indices(pixel_cluster, wild_indices[j]))[0] + start_x
                  wild_y = (array_indices(pixel_cluster, wild_indices[j]))[1] + start_y
                  num_times_flagged_as_spike[wild_x, wild_y] += 1
               endfor
            endif
         
         endfor  ; next y pixel
      endfor     ; next row

         
      ; where spikes have been identified, make image pixel invalid   
      ; Note: currently requires a pixel to be identified as a spike in all
      ; clusters checked
      spike_indices = where(num_times_flagged_as_spike eq num_times_checked, $
                            num_spikes)
   
      if num_spikes gt 0 then begin
    
         (*(images[i].image))[spike_indices] = !VALUES.F_NAN 
         (*(images[i].image_unc))[spike_indices] = !VALUES.F_NAN 
   
         if get_debug_level() ge 5 then begin
            for spike_index=0, num_spikes-1 do begin
               message, "spike found in science image with timetag " + $
                       string(images[i].image_tlm_timestamp, format='(i0)'),  " at: " + $
                       string(array_indices(*(images[i].image), spike_indices[spike_index]), format='(i4)'), /info
            endfor
         endif
      endif      
   endfor     ; next image    

   return, images
end