;+
; $Id: hi_remove_starfield.pro,v 1.1 2008/01/21 16:11:57 bewsher Exp $
;
; Project: STEREO-SECCHI
;
; Name: hi_remove_starfield
;
; Purpose: To remove stars from HI data (though it should work for
;          other types of data) by identifying strong sources in the
;          del^2 field. Only pixels identified as stars are
;          reconstructed, other pixels retain their values. Note that
;          the function returns a floating point array regardless of
;          input data type
;
; Category: STEREO, SECCHI, HI
;
; Explanation: 
;
; Syntax: recon=hi_remove_starfield(data)
;
; Example: 
;
; Inputs: data - 2d image 
;
; Opt. Inputs: thresh - this is a threshold for the strength of the del^2
;              sources. Sources stronger than this are removed. The
;              best value to use will vary depending on instrument
;              (HI1/2) and exposure time. A rule of thumb for deciding 
;              on a value for this keyword is to plot your uncorrected
;              image with your bytscaling, then set thresh to around
;              10% of your max value. For example, if you bytscale
;              your data between 0 and 1000, then try using a thresh
;              of 100.
;
;              You may decide to refine this further after seeing
;              initial results. Beware of making this value too small, 
;              as this may start to nibble away at coronal structure
;              that you may want to keep. A method to check whether
;              this is happening is to calculate the residuals
;              (data-recon) and view this as an overexposed image
;              (min=-10, max=10 say), this should indicate whether the
;              coronal structure has been affected, or whether only
;              stars (or data spikes) have been removed.  
;
;              edgefix - this keyword also has a go at cleaning up the
;                     edge pixels of the data. Spikes at the edge can
;                     lead to bleeding into the data, particularly
;                     with small data sets.
;
; Outputs: recon - reconstructed data with stars removed
;
; Opt. Outputs: None
;
; Keywords: None
;
; Calls: hi_delsq, hi_index_peakpix, hi_sor2, hi_edge_fix
;
; Common: None
;
; Restrictions: None 
;
; Side effects: None
;
; Prev. Hist.: None
;
; History: Written 27 February 2007 by Daniel Brown (DOB)
;          23-Mar-07, documentation added (DOB)
;          05-May-07, /edgefix keyword added
;
; Contact: Daniel Brown (dob@aber.ac.uk)
;
; $Log: hi_remove_starfield.pro,v $
; Revision 1.1  2008/01/21 16:11:57  bewsher
; First release of hi_remove_starfield and associated code
;
;-
function hi_remove_starfield,data,thresh=thresh,edgefix=edgefix

; set a default threshold
if (n_elements(thresh) eq 0) then thresh=1.0

; calculate the del^2 field of the data
print,'Calculating del^2 field'
divim=hi_delsq(data)

; look for pixels in the |del^2| field above the threshold
; along with their immediate neighbours
print,'Searching for out of range pixels'
ind=hi_index_peakpix(divim,thresh=thresh)

data2=data
if (keyword_set(edgefix)) then data2=hi_edge_fix(data,thresh=thresh)

if (ind(0,0) ne -999) then begin
    print,'Correcting ',strcompress(string(n_elements(ind(0,*))),/remove_all),' pixels'

; set all selected pixels in the del^2 field to 0 (this is probably
; equivalent to a 2D linear fit for reconstructing pixels. This could
; be modified to something more complex (which was the original
; intention), however, this seems to give satisfactory results
; (certainly, any alternate replacement values would have to be below
; the threshold, and ideally much smaller anyway) 
    divim(ind(0,*),ind(1,*))=0.
    recon=hi_sor2(data2,divim,ind)
endif else begin
    print,'no pixels require reconstruction'
    recon=data2
endelse

; return reconstructed image
return,recon
end