;+
; $Id: hi_edge_fix.pro,v 1.1 2008/01/21 16:14:31 bewsher Exp $
;
; Project: STEREO-SECCHI
;
; Name: hi_edge_fix
;
; Purpose: to apply a star removal algorithm to the edge of an
;          image. This uses a 1D variant of the main
;          hi_remove_starfield algorithm, but can be called by that
;          routine by a keyword argument.
;
; Category: STEREO, SECCHI, HI
;
; Explanation: If a star occurs on the edge of an image, it will not
;              be flagged for removal by the main hi_remove_starfield
;              algorithm. This is due to the nature of calculating
;              del^2. This may cause bleeding into the image during
;              reconstruction of pixels. Calling this function first
;              will minimise this effect, this can be done with a
;              keyword argument in hi_remove_starfield.
;
; Syntax: dout = hi_edge_fix(data, thresh=thresh)
;
; Example: 
;
; Inputs: data - 2d image 
;
; Opt. Inputs: thresh - a threshold for the the intensity of spike to
;                       be identified as a star. 
;
; Outputs: dout - a 2d image with cleaned edges
;
; Opt. Outputs: None
;
; Keywords: None
;
; Calls: None
;
; Common: None
;
; Restrictions: None 
;
; Side effects: None
;
; Prev. Hist.: None
;
; History: Written 16 July 2007 by Daniel Brown (DOB)
;
; Contact: Daniel Brown (dob@aber.ac.uk)
;
; $Log: hi_edge_fix.pro,v $
; Revision 1.1  2008/01/21 16:14:31  bewsher
; First release of hi_remove_starfield and associated code
;
;-
function hi_edge_fix,data,thresh=thresh

if (n_elements(thresh) eq 0) then thresh=1.0

nits=100
w=1.5

sz=size(data)
dout=data

ln = [reform(data(*,0)), reform(data(sz(1)-1,1:sz(2)-2)), reverse(reform(data(*,sz(2)-1))), reverse(reform(data(0,1:sz(2)-2)))]

nn=n_elements(ln)

dln=fltarr(nn)

dln(0) = ln(nn-1) - 2.*ln(0) + ln(1)
for i=1,nn-2 do begin
    dln(i) = ln(i-1) - 2.*ln(i) + ln(i+1)
endfor
dln(nn-1) = ln(nn-2) - 2.*ln(nn-1) + ln(0)


npix=0

tst1=abs(dln(nn-1)) gt thresh
tst2=abs(dln(0)) gt thresh
tst3=abs(dln(1)) gt thresh
if (tst1 or tst2 or tst3) then npix=npix+1

for i=1,nn-2 do begin
    tst1=abs(dln(i)) gt thresh
    tst2=abs(dln(i-1)) gt thresh
    tst3=abs(dln(i+1)) gt thresh
    if (tst1 or tst2 or tst3) then npix=npix+1
endfor

tst1=abs(dln(nn-2)) gt thresh
tst2=abs(dln(nn-1)) gt thresh
tst3=abs(dln(0)) gt thresh
if (tst1 or tst2 or tst3) then npix=npix+1


if (npix gt 0) then begin

    idx=intarr(npix)
    n=0
    
    tst1=abs(dln(nn-1)) gt thresh
    tst2=abs(dln(0)) gt thresh
    tst3=abs(dln(1)) gt thresh
    if (tst1 or tst2 or tst3) then begin
        idx(n)=0
        n=n+1
    endif

    for i=1,nn-2 do begin
        tst1=abs(dln(i)) gt thresh
        tst2=abs(dln(i-1)) gt thresh
        tst3=abs(dln(i+1)) gt thresh
        if (tst1 or tst2 or tst3) then begin
            idx(n)=i
            n=n+1
        endif
    endfor

    tst1=abs(dln(nn-2)) gt thresh
    tst2=abs(dln(nn-1)) gt thresh
    tst3=abs(dln(0)) gt thresh
    if (tst1 or tst2 or tst3) then begin
        idx(n)=nn-1
        n=n+1
    endif


    dln(idx)=0

    for j=1,nits do begin
        for k=0,npix-1 do begin

            ii1=idx(k)-1
            ii2=idx(k)+1
            if (ii1 lt 0) then ii1=nn-1
            if (ii2 gt nn-1) then ii2=0

            dum = 0.5*(ln(ii1) + ln(ii2) - dln(idx(k)))
            dum2 = w*dum + (1.0-w)*ln(idx(k))
            ln(idx(k)) = dum2

        endfor
    endfor


    i1=sz(1)
    i2=sz(1) + sz(2)-2
    i3=2*sz(1) + sz(2)-2
    i4=2*sz(1) + 2*sz(2)-4

    dout(*,0) = ln(0:i1-1)
    dout(sz(1)-1,1:sz(2)-2) = ln(i1:i2-1)
    dout(*,sz(2)-1) = reverse(ln(i2:i3-1))
    dout(0,1:sz(2)-2) = reverse(ln(i3:i4-1))

endif

return,dout
end