function gauss_remove, img, wsize = wsize, threshold = threshold, show $
  = show, order = order, residual = residual
;+
; NAME:
;	GAUSS_REMOVE
;
;
; PURPOSE:
;	Clean up an image by subtracting discrete gauassian features
;	from it.
;
;
; CATEGORY:
;	Utils
;
;
; CALLING SEQUENCE:
;	imgn = gauss_remove(img[, wsize=wsize, threshold=threshold, $
;		show=show, order=order, residual=residual])
;
;
; INPUTS:
;	img	float	The image to be cleaned, this must be a named
;			variable as ithe cleaing is done in place.
;
;
; KEYWORD PARAMETERS:
;	wsize	int	The size of the window to use for the fitting
;			(default 25x25)
;	threshold float The theshold for subtracting a feature.
;	show	int	If present and non-zero display extra
;			information about the fitting, if explicitly
;			zero, then suppress progress info.
;	order	int	The order of the polynomial fit to precede the
;			gaussian fit (default 2)
;	residual float	A named variable to hold the blemishes subtracted.
;
; OUTPUTS:
;	imgn	float	The image with the "blemishes" removed
;
;
; MODIFICATION HISTORY:
;	Original: 28/7/04; SJT
;	Converted to function and renamed: 30/7/04; SJT
;-


if not keyword_set(wsize) then wsize = [25, 25]
if n_elements(wsize) eq 1 then wsize = replicate(wsize, 2)

if not keyword_set(threshold) then threshold = 2.
if n_elements(show) ne 0 and not keyword_set(show) then noplot = 1b $
else noplot = 0b
if n_elements(order) eq 0 then order = 2

steps = floor(wsize/3) > 1
wb = floor(wsize/4)
wt = wsize-wb

sz = size(img)
residual =  make_array(size = sz)
imgn = img

if keyword_set(show) then window, 3
for iy = steps[1], sz[2]-steps[1], steps[1] do begin
    iy0 = (iy-steps[1]) > 0
    iy1 = (iy+steps[1]) < (sz[2]-1)
    if iy1-iy0 lt 3 then continue
    for ix = steps[0], sz[1]-steps[0], steps[0] do begin
        ix0 = (ix-steps[0]) > 0
        ix1 = (ix+steps[0]) < (sz[1]-1)
        if ix1-ix0 lt 3 then continue
        subi = imgn[ix0:ix1, iy0:iy1]
        locs = where(finite(subi), comp = nlocs, nf, ncomp = nn)
        if nf lt n_elements(subi)/2. then continue
        avg = mean(subi, /nan)
        mx = max(subi, /nan, min = mn)
        if mx - mn lt 2*threshold then continue
        if nn ne 0 then subi[nlocs] = avg
        sf = sfit(subi, order)
        subi -= sf
        avg1 = mean(subi, /nan)
        mx1 = max(subi, /nan, min = mn1)
        if mx1 - avg1 gt avg1 - mn1 then ft = gauss2dfit(subi, a, /tilt) $
        else  ft = gauss2dfit(subi, a, /tilt, /negative)
        if (a[4] lt wb[0] or a[4] gt wt[0] and $ ; Only subtract peaks
            a[5] lt wb[1] or a[5] gt wt[1]) then continue ; near the
                                ; middle of the box
        if a[2]*a[3] gt wsize[0]*wsize[1]/4. and (a[2] < a[3]) gt $
          min(wsize)/5. then continue ; Too big a structure
        if keyword_set(show) then begin
            surface, ft
            empty
        endif
        ft -= median(ft)
        if nn ne 0 then ft[nlocs] = 0
        imgn[ix0:ix1, iy0:iy1] -= ft
        residual[ix0:ix1, iy0:iy1] += ft
        if not noplot then print, ix, iy,  nf, nn
    endfor
endfor
if not noplot then begin
    window, 0, xs = sz[1], ys = sz[2]
    tv, bytscl(img, min = -1, max = 1)
    window, 1, xs = sz[1], ys = sz[2]
    tv, bytscl(imgn, min = -1, max = 1)
    window, 2, xs = sz[1], ys = sz[2]
    tv, bytscl(residual, min = -1, max = 1)
endif

return, imgn

end