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.
;
; Modified 25 Sept 1992 by DPS
; Each pixel indexed by 'toobig' is set to the median of the 3x3
; neighbourhood centered on the offending pixel.  Pixels on the boundaries
; of the image are not changed, since the neighbourhood is undefined.
;
l=256.
col=REBIN(LINDGEN(l),l,l)
row=TRANSPOSE(col)
w=where(SQRT((col-131)^2. + (row-125)^2.) LT 146)
imgr=intarr(256,256)
imgr(w)=image(w)
dim=SIZE(imgr)
g=ABS(imgr-SHIFT(imgr,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(imgr-SHIFT(imgr,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(imgr-SHIFT(imgr,-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(imgr-SHIFT(imgr,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)
r=toobig/dim(1)
c=toobig MOD dim(1)
tmp=WHERE(((dr-1 LT c) AND (c LT dim(1)-dr)) AND $
          ((dr-1 LT r) AND (r LT dim(2)-dr)))
IF tmp(0) NE -1 THEN BEGIN
    toobig=toobig(tmp)
    r=r(tmp) & c=c(tmp)
    PRINT,'clean :',STRTRIM(N_ELEMENTS(toobig)),' values changed'
    changed=imgr	; copy input image
    FOR i=0,N_ELEMENTS(toobig)-1 DO $
        changed(toobig(i))=MEDIAN(imgr(c(i)-1:c(i)+1,r(i)-1:r(i)+1))
    RETURN,changed
ENDIF ELSE BEGIN
    PRINT,'clean : 0 values changed'
    RETURN,imgr
ENDELSE
END