FUNCTION clean,image,dr,ns,ng
;
; This procedure is intended to clean a CCD image of random high
; counts, or "cosmic ray" events. It is based on an algorithm from
; "Preprocessing of CCD images", by Peter Linde of the Lund Observatory,
; Lund, Sweden.
; 
; The algorithm calculates image gradients over a neighbourhood size
; dr, and zeroes the vectors where image wraparound has produced
; meaningless gradients. A gradient threshold for each gradient image,
; equal to the mean gradient plus ns times the standard deviation in
; the gradient, is calculated for each gradient image. The maximum
; of these four gradients is taken as the gradient threshold for
; identifying abnormal gradients. The number of gradients (out of
; 4) at each pixel that exceed this threshold is determined, and
; those pixels that have more than ng high gradients are set to the
; median pixel value of the input image. The modified image is
; returned.
;
dim=SIZE(image)
level=MEDIAN(image)		; value to substitute for high counts
g=ABS(image-SHIFT(image,dr,0))	; calculate gradient to right and
g(dr-1,*)=0			; zero meaningless gradients
sd=STDEV(g,m)&gt1=m+ns*sd	; calculate gradient threshold
g1big=WHERE(g GT gt1)		; find big gradients.
IF g1big(0) NE -1 THEN big1=g(g1big) ELSE big1=-1
; if some big gradients found, big1 is a list of them, otherwise empty
g=ABS(image-SHIFT(image,0,dr))	; repeat process for gradient upward ...
g(*,dim(2)-dr)=0
sd=STDEV(g,m)&gt2=m+ns*sd
g2big=WHERE(g GT gt2)
IF g2big(0) NE -1 THEN big2=g(g2big) ELSE big2=-1
g=ABS(image-SHIFT(image,-dr,0))	; repeat process for gradient to left ...
g(dim(1)-dr,*)=0
sd=STDEV(g,m)&gt3=m+ns*sd
g3big=WHERE(g GT gt3)
IF g3big(0) NE -1 THEN big3=g(g3big) ELSE big3=-1
g=ABS(image-SHIFT(image,0,-dr))	; repeat process for gradient downward.
g(*,dr-1)=0
sd=STDEV(g,m)&gt4=m+ns*sd
g4big=WHERE(g GT gt4)
IF g4big(0) NE -1 THEN big4=g(g4big) ELSE big4=-1
gradthresh=MAX([gt1,gt2,gt3,gt4])	; select the largest gradient
vbig=WHERE(big1 GT gradthresh)		; restrict big gradients to
IF vbig(0) NE -1 THEN g1big=g1big(vbig) ELSE g1big=-1; those greater than gradthresh
vbig=WHERE(big2 GT gradthresh)
IF vbig(0) NE -1 THEN g2big=g2big(vbig) ELSE g2big=-1
vbig=WHERE(big3 GT gradthresh)
IF vbig(0) NE -1 THEN g3big=g3big(vbig) ELSE g3big=-1
vbig=WHERE(big4 GT gradthresh)
IF vbig(0) NE -1 THEN g4big=g4big(vbig) ELSE g4big=-1
nbig=BYTARR(dim(1),dim(2))	; no. of big gradients at each pixel
IF g1big(0) NE -1 THEN nbig(g1big)=nbig(g1big)+1
IF g2big(0) NE -1 THEN nbig(g2big)=nbig(g2big)+1
IF g3big(0) NE -1 THEN nbig(g3big)=nbig(g3big)+1
IF g4big(0) NE -1 THEN nbig(g4big)=nbig(g4big)+1
toobig=WHERE(nbig GE ng)&change=(toobig(0) NE -1)
IF change THEN BEGIN
    PRINT,'clean :',STRCOMPRESS(STRING(N_ELEMENTS(toobig))),' values changed'
    changed=image&changed(toobig)=level	; copy input image, correct high-count pixels
    RETURN,changed
ENDIF ELSE BEGIN
    PRINT,'clean : 0 values changed'
    RETURN,image
ENDELSE
END