;+
; $Id: hi_delsq.pro,v 1.2 2009/06/18 13:54:49 crothers Exp $
;
; Project: STEREO-SECCHI
;
; Name: hi_delsq
;
; Purpose: To calculate the del^2 field of an image
;
; Category: STEREO, SECCHI, HI
;
; Explanation: 
;
; Syntax: del2=hi_delsq(data)
;
; Example: 
;
; Inputs: data - 2d image 
;
; Opt. Inputs: LEAST - least number of valid neighbours [default 4]
;
; Outputs: del2 - del^2 field
;
; Opt. Outputs: None
;
; Keywords: None
;
; Calls: None
;
; 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)
;
; Contact: Daniel Brown (dob@aber.ac.uk)
;
; $Log: hi_delsq.pro,v 
; Revision 1.1  2008/01/21 16:13:09  bewshe
; First release of hi_remove_starfield and associated cod
;
;-
function hi_delsq,dat,least=least

if n_elements(least) eq 0 then least=4

; find extent of data
sz=size(dat)
nx=sz(1)
ny=sz(2)

; create a working array
dhj=dat*0.0
ctr=intarr(nx,ny)

; forward addressing

d0=dat[0:nx-2,*]-dat[1:nx-1,*]

map=finite(d0)
index=where(~map,nindex)
if nindex gt 0 then d0[index]=0

dhj[0:nx-2,*]=dhj[0:nx-2,*]+d0
ctr[0:nx-2,*]=ctr[0:nx-2,*]+map

; backward addressing

d0=dat[1:nx-1,*]-dat[0:nx-2,*]

map=finite(d0)
index=where(~map,nindex)
if nindex gt 0 then d0[index]=0

dhj[1:nx-1,*]=dhj[1:nx-1,*]+d0
ctr[1:nx-1,*]=ctr[1:nx-1,*]+map

; downard addressing

d0=dat[*,0:nx-2]-dat[*,1:nx-1]

map=finite(d0)
index=where(~map,nindex)
if nindex gt 0 then d0[index]=0

dhj[*,0:nx-2]=dhj[*,0:nx-2]+d0
ctr[*,0:nx-2]=ctr[*,0:nx-2]+map

 ; upward addressing

d0=dat[*,1:nx-1]-dat[*,0:nx-2]

map=finite(d0)
index=where(~map,nindex)
if nindex gt 0 then d0[index]=0

dhj[*,1:nx-1]=dhj[*,1:nx-1]+d0
ctr[*,1:nx-1]=ctr[*,1:nx-1]+map

index=where(ctr lt least,nindex)
if nindex gt 0 then dhj[index]=!values.f_nan

; return field
return,dhj

end