PRO CDS_CLEAN_IMAGE, IMAGE, MISSING=MISSING, NOMISSING=NOMISSING, $
		MINSIG=MINSIG
;+
; Project     :	SOHO - CDS
;
; Name        :	CDS_CLEAN_IMAGE
;
; Purpose     :	Clean cosmic rays & missing pixels from NIS raster images
;
; Category    :	Class3, Image-display
;
; Explanation :	This routine removes cosmic rays and missing pixels from NIS
;		raster images.
;
;		Cosmic rays are recognized as pixels which are very high
;		compared to neighboring exposures.  Thus, this procedure can
;		only be used on images which are built up from rastering the
;		slit across the sun.  The first dimension must represent the
;		raster dimension, and the second dimension must represent
;		points along the slit.  Spectra are not supported.
;
;		The algorithm used to recognize cosmic rays is far from
;		perfect.  It is recommended that this routine be used only for
;		image display purposes, and not for analysis.
;
; Syntax      :	CDS_CLEAN_IMAGE, IMAGE  [, MISSING=MISSING ]
;
; Examples    :	A = READCDSFITS('s2707r00.fits')
;		IMAGE = AVERAGE(DETDATA(A,1),1,MISSING=-1)
;		SETFLAG, MISSING=-1
;		CDS_CLEAN_IMAGE, IMAGE
;		EXPTV, IMAGE
;
; Inputs      :	IMAGE	= The image to be cleaned.
;
; Opt. Inputs :	None.
;
; Outputs     :	IMAGE	= The cleaned image is returned in the same parameter.
;
; Opt. Outputs:	None.
;
; Keywords    :	MISSING	= Value representing missing pixels.  If not passed,
;			  then the value set by SETFLAG is used, or the routine
;			  assigns it's own value, to use for filtering out
;			  pixels containing cosmic ray hits.
;
;		NOMISSING = If set, then only clean cosmic rays--don't try to
;			    clean up missing pixels.  However, the MISSING
;			    value can still be passed, through either the
;			    keyword or SETFLAG, to mark pixels that really are
;			    missing.
;
;		MINSIG	= Minimum value to use when estimating the standard
;			  deviation in the pixels.
;
; Calls       :	GET_IM_KEYWORD, FILL_MISSING
;
; Common      :	None.
;
; Restrictions:	None, except those noted above.
;
; Side effects:	None.
;
; Prev. Hist. :	None.
;
; History     :	Version 1, 17-May-1996, William Thompson, GSFC
;		Version 2, 05-Jun-1996, William Thompson, GSFC
;			Make more robust when roundoff errors occur
;               Version 3, 20-May-1997, Zarro, GSFC, stopped abrupt
;                       return caused by MESSAGE
;		Version 4, 14-Aug-1997, William Thompson, GSFC
;			Added keyword NOMISSING
;		Version 5, 17-Apr-2000, William Thompson, GSFC
;			Extended version 3 change to other MESSAGE calls
;
; Contact     :	WTHOMPSON
;-
;
	ON_ERROR, 2
;
;  Check the input parameters.
;
	IF N_PARAMS() NE 1 THEN BEGIN
	    MESSAGE, 'Syntax:  CDS_CLEAN_IMAGE, IMAGE  [, MISSING=MISSING]', $
		    /CONTINUE
	    RETURN
	ENDIF
;
;  Make sure that the input image is two-dimensional.
;
	SZ = SIZE(IMAGE)
	IF SZ(0) NE 2 THEN BEGIN
	    MESSAGE, 'IMAGE must have two dimensions', /CONTINUE
	    RETURN
	ENDIF
	NX = SZ(1)
	NY = SZ(2)
	IF (NX LT 3) OR (NY LT 3) THEN BEGIN
	    MESSAGE, 'IMAGE must be at least 3x3 in size', /CONTINUE
	    RETURN
        ENDIF
;
;  Get the missing pixel flag value.
;
	GET_IM_KEYWORD, MISSING, !IMAGE.MISSING
;
;  If the MISSING keyword wasn't set then use the smallest value in the array,
;  minus 1.  Otherwise, store the current value for later use.
;
	IF N_ELEMENTS(MISSING) NE 1 THEN MISSING = MIN(IMAGE) - 1 ELSE	$
		KMISSING = MISSING
;
;  If the NOMISSING keyword was set, then use the default value anyway.
;
	TWO_MISSINGS = KEYWORD_SET(NOMISSING) AND (N_ELEMENTS(KMISSING) EQ 1)
	IF TWO_MISSINGS THEN MISSING = MIN(IMAGE) - 1
;
;  Make sure there are some good pixels.
;
	IF TWO_MISSINGS THEN GOOD = IMAGE NE KMISSING ELSE	$
		GOOD = IMAGE NE MISSING
	IF TOTAL(GOOD) EQ 0 THEN BEGIN
	    MESSAGE, 'No pixels present', /CONTINUE
	    RETURN
	ENDIF
;
;  Form the boxcar total and total square of the image over three vertical
;  pixels at each point.
;
	TEMP = 1.*IMAGE
	BAD = 1B - GOOD
	W_BAD = WHERE(BAD, N_BAD)
	IF N_BAD GT 0 THEN TEMP(W_BAD) = 0
	TOT = TEMP
	ASQ = TEMP^2
	NAV = GOOD
	TOT(*,0:NY-2) = TOT(*,0:NY-2) + TEMP(*,1:*)
	ASQ(*,0:NY-2) = ASQ(*,0:NY-2) + TEMP(*,1:*)^2
	NAV(*,0:NY-2) = NAV(*,0:NY-2) + GOOD(*,1:*)
	TOT(*,1:*) = TOT(*,1:*) + TEMP(*,0:NY-2)
	ASQ(*,1:*) = ASQ(*,1:*) + TEMP(*,0:NY-2)^2
	NAV(*,1:*) = NAV(*,1:*) + GOOD(*,0:NY-2)
;
;  Step through each column in the image, and see if it's more than 3 sigma
;  above its (wrapped around) neighbors.
;
	TAV = TOT
	TAV(1:NX-2,*) = (TOT(0:NX-3,*) + TOT(2:*,*)) /	$
		((NAV(0:NX-3,*) + NAV(2:*,*)) > 1)
	TAV(0,*)    = (TOT(NX-1,*) + TOT(1,*)) / ((NAV(NX-1,*) + NAV(1,*)) > 1)
	TAV(NX-1,*) = (TOT(NX-2,*) + TOT(0,*)) / ((NAV(NX-2,*) + NAV(0,*)) > 1)
;
	TSQ = ASQ
	TSQ(1:NX-2,*) = (ASQ(0:NX-3,*) + ASQ(2:*,*)) /	$
		((NAV(0:NX-3,*) + NAV(2:*,*)) > 1)
	TSQ(0,*)    = (ASQ(NX-1,*) + ASQ(1,*)) / ((NAV(NX-1,*) + NAV(1,*)) > 1)
	TSQ(NX-1,*) = (ASQ(NX-2,*) + ASQ(0,*)) / ((NAV(NX-2,*) + NAV(0,*)) > 1)
;
	SIG = SQRT((TSQ - TAV^2) > 0)
	IF N_ELEMENTS(MINSIG) EQ 1 THEN SIG = SIG > MINSIG
	W = WHERE((TEMP GT TAV+3*SIG) AND (SIG GT 0), N_CR)
	IF N_CR GT 0 THEN BAD(W) = 1
;
;  If the NOMISSING keyword was passed, then save the positions of all pixels
;  that really are missing.
;
	IF TWO_MISSINGS THEN W_MISSING = WHERE(IMAGE EQ KMISSING, C_MISSING)
;
;  Set all the cosmic rays found to the missing pixel flag value.
;
	W_BAD = WHERE(BAD, N_BAD)
	IF N_BAD GT 0 THEN BEGIN
		IMAGE(W_BAD) = MISSING
;
;  Step through the rows and replace any missing pixels.
;
		FOR IY = 0,NY-1 DO BEGIN
			TEMP = IMAGE(*,IY)
			FILL_MISSING, TEMP, MISSING
			IMAGE(*,IY) = TEMP
		ENDFOR
	ENDIF
;
;  If the NOMISSING keyword was passed, then reset the missing pixels.
;
	IF TWO_MISSINGS THEN IF C_MISSING GT 0 THEN IMAGE(W_MISSING) = KMISSING
;
	RETURN
	END