PRO CLEAN_SPIKE, ARR_IN, ARR_OUT=ARR_OUT, SPIKES=SPIKES, CUTOFF=CUTOFF, $
	BCKGRND=BCKGRND, MISSING=MISSING, REPLACE=REPLACE, $
	PRETTY=PRETTY, FILL=FILL, NOINFO=NOINFO
;+
; Project     : SOHO - CDS
;
; Name        : CLEAN_SPIKE
;
; Purpose     : Removes cosmic rays from a 3-D CDS array
;
; Category    : General utility.
;
; Explanation : Takes a (LAMB,X,Y) CDS array and looks in the X-direction for 
;		"spikes" (ie, pixel is greater by a factor CUTOFF than it's 
;		two X-neighbours, and must have an intensity greater than 
;		BCKGRND). For confirmation of the spike identification, the 
;		program does the same checks in the LAMB and Y directions. If 
;		at least one of these confirms the spike, then it is classed 
;		as a cosmic ray. The cosmic ray can be set to 'MISSING' 
;		(useful for analysis) or to the average of the two 
;		neighbouring X-pixels (useful for pretty pictures).
;
;		Some cosmic rays have a 2-pixel spike in the Y direction and 
;		so a check is made in a 4th pixel adjacent to one of the 
;		neighbouring pixels to account for this situation.
;
;		The BCKGRND value is required to prevent noise spikes in the 
;		background from being removed.
;
; Syntax      : CLEAN_SPIKE, ARR_IN, [ ARR_OUT=ARR_OUT, SPIKES=SPIKES, 
;			CUTOFF=CUTOFF, BCKGRND=BCKGRND, MISSING=MISSING,
;			/REPLACE, /PRETTY ]
;
; Examples    : CLEAN_SPIKE, ARR_IN, /REPLACE
;
;		CLEAN_SPIKE, ARR_IN, ARR_OUT=ARR_OUT, SPIKES=SPIKES
;
; Inputs      : ARR_IN	= A 3-D data array, e.g., a data window extracted 
;			  using GT_WINDATA
;
; Opt. Inputs :	CUTOFF	= Number for determining definition of spike. 
;
;		BCKGRND	= Number specifying level of background.
;
;		MISSING = The value of possible missing data (like
;                         already-identified cosmic rays), and the value used
;                         to flag the cosmic ray pixels. If this value is not
;                         supplied, the pixel value is replaced by the average
;			  of the two neighbouring X pixels.
;
; Outputs     : None.
;
; Opt. Outputs: ARR_OUT	= The cleaned data array.
;
; 		SPIKES	= An array of the same dimension of ARR_IN containing 
;			  the locations of all the cosmic rays. A cr is 
;			  denoted by a '1'.
;
; Keywords    : REPLACE = If set, this keyword replaces ARR_IN with the 
;			  cleaned array.
;
;		PRETTY	= Sets more strict parameters for determining cosmic 
;			  rays, and hence makes the images prettier!
;
;		FILL    = Fills in the cosmic ray with the average of the 
;			  two neighbouring X pixels. Over-rides the keyword
;			  MISSING if this is set.
;
;		NOINFO	= Suppresses information about number of identified
;			  cosmic ray pixels
;
; Calls       : Y_CHECKER, RAY_NEIGHBOUR, FMEDIAN
;
; Common      : None.
;
; Restrictions: ARR_IN must be 3-D and each dimension 
;		must be greater than 2. Will not work on 2-D images!
;
;		Data needs to be de-biased first.
;
;		For some large cosmic rays two runs of CLEAN_SPIKE may be 
;		required to remove the cosmic ray.
;
; Side effects: None.
;
; Prev. Hist. : None.
;
; History     : Ver. 1 - PRY, 19/9/96
;		Ver. 2 - PRY, 4/12/96 have added an estimation of the local 
;				background using the FMEDIAN routine.
;		Ver. 3 - PRY, 24/9/97 corrected problem with keyword FILL
;               Ver. 4 - PRY, 20/11/97 routine displays list of parameters if 
;			        none are given
;               Ver. 5 - PRY, 20/11/97 corrected bug in the use of the
;		     		REPLACE keyword
;
; Contact     :  Peter Young,  Cambridge University (pry10@damtp.cam.ac.uk)
;-
;

IF N_PARAMS() LT 1 THEN BEGIN
  PRINT,'Use: IDL> clean_spike, arr_in, [ arr_out=arr_out, spikes=spikes, 
  PRINT,'                          cutoff=cutoff, bckgrnd=bckgrnd, 
  PRINT,'                          missing=missing, /replace, 
  PRINT,'                          /pretty, /fill, /noinfo ]
  RETURN
ENDIF

; ---------------
; Set the parameters that determine a cosmic ray
;
IF N_ELEMENTS(cutoff) eq 0 THEN cutoff=.6
IF N_ELEMENTS(bckgrnd) eq 0 THEN bckgrnd=8
;
IF KEYWORD_SET(pretty) THEN BEGIN
  cutoff=.4
  bckgrnd=5
ENDIF
; ---------------

arb_amount=0.1

dim=size(arr_in) & ARR_DESP=arr_in-arr_in & SPIKES=ARR_DESP

any_dropouts=0
IF TOTAL(ARR_in EQ -100) GT 0 THEN any_dropouts=1

l_siz=dim(1)-1
x_siz=dim(2)-1
y_siz=dim(3)-1


;
; There are three options for the replacement value of the cosmic ray,
; depending on the keywords used:
;     1. MISSING set and FILL not set -- use MISSING value
;     2. FILL set  -- use average of neighbouring x-pixels
;     3. MISSING not set and FILL not set  -- replace with a value of 100
; These three options will be noted by the value of replace.
;
IF (N_ELEMENTS(missing) NE 0) AND (NOT KEYWORD_SET(fill)) THEN rep_CR=1
IF KEYWORD_SET(fill) THEN rep_CR=2
IF (N_ELEMENTS(missing) EQ 0) AND (NOT KEYWORD_SET(fill)) THEN BEGIN
  rep_CR=3
  PRINT,'Spikes will be replaced with a value of -100. To replace spikes'
  PRINT,'with an average of neighbouring pixels, use the /FILL keyword.'
  PRINT,''
ENDIF

FOR p=0,1 DO BEGIN
  FOR i=0,l_siz DO BEGIN
    lp=p * l_siz + i * (1 - 2*p)
    img=REFORM(ARR_in(lp,*,*))
    med_img=fmedian(img,5,5)
    rays=img-img   ; a blank array of the same size as img

    FOR q=0,1 DO BEGIN
      FOR j=0,y_siz DO BEGIN
        yp=q * y_siz + j * (1 - 2*q)
        ;
        row=REFORM(img(*,yp))
        med_row=REFORM(med_img(*,yp))
	;
        ray_neighbour,row,av,arb_amount,any_dropouts
        ind=WHERE( ( (row/av - 1) GT cutoff ) AND ((row-med_row) GT bckgrnd) )
        IF ind(0) NE -1 THEN BEGIN
          n=N_ELEMENTS(ind)
          rem_ind=-1
          FOR k=1,n DO BEGIN
            index=ind(k-1)
;
; First check in the y and lambda directions that the feature is a cosmic
; ray. Will apply the same cutoff as for the x direction.
;
; y_ray and l_ray are logical variables, set to true if cosmic ray is there.
;
            l_ray=0
;
; do lambda first
; ---------------
            CASE lp OF
;
            0: BEGIN
                 l_col=REFORM(ARR_IN(lp:lp+1,index,yp))
		 ray_neighbour,l_col,l_col_av,arb_amount,any_dropouts
	         IF ( l_col(1)/l_col_av(1) - 1) GT cutoff THEN l_ray=1
               END

            l_siz: BEGIN
	             l_col=REFORM(ARR_IN(lp-1:lp,index,yp))
		     ray_neighbour,l_col,l_col_av,arb_amount,any_dropouts
	             if ( l_col(1)/l_col_av(1) - 1) gt cutoff THEN l_ray=1
                   END

            else: BEGIN
                    l_col=REFORM(ARR_IN(lp-1:lp+1,index,yp))
	            ray_neighbour,l_col,l_col_av,arb_amount,any_dropouts
	            IF ( l_col(1)/l_col_av(1) - 1) GT cutoff THEN l_ray=1
                  END
;
            ENDCASE
;
; do y
; ----
            if l_ray eq 0 THEN BEGIN
	      col=REFORM(img(index,*))
              y_checker,col,yp,y_siz,cutoff,y_ray,arb_amount,any_dropouts
            ENDIF else y_ray=0
;
;
; If neither l_ray or y_ray have been set, THEN we reject the cosmic ray
; removal, and take index out of ind.
;
            IF (rem_ind(0) NE -1) THEN BEGIN
              IF (l_ray EQ 1) OR (y_ray EQ 1) THEN rem_ind=[rem_ind,index]
            ENDIF

            IF (rem_ind(0) EQ -1) THEN BEGIN
              IF (l_ray eq 1) OR (y_ray eq 1) THEN rem_ind=index
            ENDIF
   
          ENDFOR	; -- end of index for loop

          ind=rem_ind

          IF ind(0) NE -1 THEN BEGIN
            ;
        ;; There are 3 options as to what the cosmic ray value will be
        ;; replaced with, depending on the keywords
        ;;
            CASE rep_CR OF
           ;
            2: row(ind)=av(ind)
           ;
            1: row(ind)=missing
           ;
            3: row(ind)=-100.
           ;
            ENDCASE
            ;
            img(*,yp)=row
            rays(ind,yp)=1
          ENDIF

        ENDIF	; -- end of "is ind empty?"

      ENDFOR	; -- end of j loop

    ENDFOR	; -- end of q loop

    SPIKES(lp,*,*)=SPIKES(lp,*,*)+rays
    ARR_DESP(lp,*,*)=img

  ENDFOR	; -- end of i loop

ENDFOR		; -- end of p loop

arr_out=arr_desp
IF KEYWORD_SET(replace) THEN arr_in=arr_desp

IF KEYWORD_SET(noinfo) EQ 0 THEN BEGIN
  ind=WHERE(spikes GT 0)
  PRINT,'                        cosmic ray pixels:',N_ELEMENTS(ind)
ENDIF

END