function d2thresh, image, threshold = threshold, median = median, $
                   kernel = kernel, nokill = nokill

;+
; D2THRESH
;	A filter based on thresholding del^2 of an image, to remove
;	spikes.
;
; Usage:
;	fimg = d2thresh(image)
;
; Returns:
;	A filtered image
;
; Argument:
;	image	float	input	The image to filter
;
; Keywords:
;	threshold float	input	The level of del^2(image) to consider
;				bad (default=3)
;	median	int	input	The size of median to use for
;				replacing the flagged values
;	kernel	int	input	Either a kernel for the DILATE
;				function or a code for a standard kernel
;	/nokill		input	If set, then don't kill the flagged
;				points before generating the median.
;
; History:
;	Original (inspired by the so-called star removal in the HI
;	software): 16/12/10; SJT
;-

case n_elements(kernel) of
    0: knl = replicate(1, 3, 3) 
    1: begin
        case kernel of
            0: knl = replicate(1, 3, 3)
            1: knl = [[0, 1, 0], $
                      [1, 1, 1], $
                      [0, 1, 0]]
            2: knl = [[0, 0, 1, 0, 0], $
                      [0, 1, 1, 1, 0], $
                      [1, 1, 1, 1, 1], $
                      [0, 1, 1, 1, 0], $
                      [0, 0, 1, 0, 0]]
            else: begin
                print, "Warning invalid kernel setting"
                knl = replicate(1, 3, 3)
            end
        endcase
    end
    else: knl = kernel
endcase

if ~keyword_set(threshold) then threshold = 3.0
if ~keyword_set(median) then median = 7

dimg = del_sq(image)
locs = where(dilate(abs(dimg) gt threshold, knl), nfix)
if nfix eq 0 then return, image

imgc = image
if ~keyword_set(nokill) then imgc[locs] = !values.f_nan
imgm = median(imgc, median)
imgc[locs] = imgm[locs]

return, imgc

end