function mcleaner, img, thresh, max_scale = max_scale, iscorr = $
                   iscorr, resid = resid, show = show
;+
; NAME:
;	MCLEANER
;
;
; PURPOSE:
;	Remove blemishes from images by successive median filters.
;
;
; CATEGORY:
;	Utils
;
;
; CALLING SEQUENCE:
;	cimg = mcleaner(img[, thresh])
;
;
; INPUTS:
;	img	float	The images to be cleaned.
;
;
; OPTIONAL INPUTS:
;	thresh	float	The threshold at which to consider a feature
;			to be a blemish (default 0.5)
;
;
; KEYWORD PARAMETERS:
;	max_scale	int	The largest median scale to use
;				(default 21)
;	iscorr		byte	A names variable to hold a map of
;				which pixels have been corrected.
;	resid		float	A named variable to hold the
;				corrections applied.
;	/show			If set, then show the progress of the
;				cleaning. (Assumes a window big enough
;				for 2 copies of the imput image is
;				available).
;
;
; OUTPUTS:
;	cimg	float	The cleaned image.
;
;
;
; PROCEDURE:
;	Median filters of increasing sizes from 3 up to max_scale are
;	applied. If a the difference between the original image and the
;	median-fltered version divided by 1 plus the sum of the
;	orignal and the filtered exceeds the threshold, then the
;	pixel is replaced by the median filtered version.
;
;
; MODIFICATION HISTORY:
;	Original: 9/8/04
;	Increase spacing of scales from 2 to 4: 24/11/04; SJT
;-

if not keyword_set(max_scale) then max_scale = 21

if n_params() ne 2 then thresh = .5
ximg = img
iscorr = bytarr(size(img, /dimension))
resid = fltarr(size(img, /dimension))

for jm = 3, max_scale, 4 do begin
    mimg = median(ximg, jm)
    kimg = abs(mimg-ximg)/(abs(mimg+ximg)+1)
    locs = where(kimg gt thresh, nb)
    if nb ne 0 then begin
        resid[locs] += (ximg[locs] - mimg[locs])
        ximg[locs] = mimg[locs]
        iscorr[locs] = 1b
    endif
    if keyword_set(show) then begin
        tvscl, kimg gt thresh, 1
        tv, bytscl(ximg, min = -2, max = 2), 0
        xyouts, .02, .02, /norm, string(jm, nb, format = "(I3,I5,'Pixels')")
        empty
    endif
endfor
if keyword_set(show) then begin
    wait, 1
    tvscl, iscorr, 1
    tv, bytscl(resid, min = -2, max = 2), 0
    wait, 1
endif

return, ximg

end