PRO eis_prep_clean_crs, wd, err, missing=missing, swtch=swtch, verbose=verbose, $
                      nofill=nofill, swap_xy=swap_xy

;+
; NAME:
;     EIS_PREP_CLEAN_CRS
;
; PURPOSE:
;     Flags and (optionally) replaces spikes in the EIS data array. An
;     EIS Software Note gives more details about the method.
;
; CATEGORY:
;     Hinode/EIS; calibration; cosmic rays.
;
; CALLING SEQUENCE:
;     EIS_PREP_CLEAN_CRS, Wd, Err
;
; INPUTS:
;     Wd:    A 3D intensity array containing an EIS data window. The dimensions
;            should be WAVELENGTH*Y*EXPOSURE (exposure being either
;            repeated exposures in sit-and-stare data or slit
;            positions in a raster scan). If residual warm pixels are
;            found then the WD values at their positions will be set
;            to interpolated values based on neighboring pixels. If
;            /nofill is set then their values are set to missing.
;     Err:   An error array of same size as WD. Missing data should be
;            flagged with the value MISSING. If residual warm pixels are
;            found then the value of ERR at their positions will be
;            set to missing.
;
; OPTIONAL INPUTS:
;     Missing: A value used for missing data. If not set then -100 is
;              assumed. 
;	
; KEYWORD PARAMETERS:
;     NOFILL:  If set, then the WD values for the warm pixels will be
;              set to MISSING, rather than interpolated values.
;     SWAP_XY: If set, then it is assumed that the WD array has
;              dimensions WAVELENGTH*EXPOSURE*Y and so the routine
;              will swap the exposure and y dimensions.
;     VERBOSE: If set, then information will be printed to the
;              screen. 
;
; OUTPUTS:
;     If residual warm pixels are found then the input arrays WD and
;     ERR will be modified (see above for more details).
;
; OPTIONAL OUTPUTS:
;     Swtch:   An byte array of same size as WD containing 0's and
;              1's. A 1 indicates the pixel has been flagged as
;              a residual warm pixel.
;
; CALLS:
;     EIS_PREP_MED_FILL, EIS_CR_CHECKER
;
; EXAMPLE:
;     This routine is called from within EIS_PREP and
;     EIS_CLEAN_WINDATA. Please see those routines for examples on how
;     it is used.
;
; MODIFICATION HISTORY:
;     Ver.1, 5-Nov-2019, Peter Young
;     Ver.2, 13-Nov-2019, Peter Young
;       Moved the checking to EIS_CR_CHECKER.
;-


IF n_params() LT 2 THEN BEGIN
  print,'Use:  IDL> eis_prep_clean_crs, wd, err [, swtch=, /nofill, /verbose, /swap_xy, missing= ]'
  return
ENDIF 


wdout=wd
errout=err

IF n_elements(missing) EQ 0 THEN missing=-100.


;
; This routine is intended to work with the 3D arrays sent by
; eis_prep, which are WVL*Y*X. If you send WVL*X*Y arrays (e.g., from
; eis_getwindata), then you need to use /SWAP_XY.
;
IF keyword_set(swap_xy) THEN BEGIN
  wdout=transpose(temporary(wdout),[0,2,1])
  errout=transpose(temporary(errout),[0,2,1])
ENDIF 

s=size(wdout)
IF s[0] EQ 2 THEN nexp=1 ELSE nexp=s[3]
;
s=size(wdout,/dim)
swtch=make_array(s,/byte)

nlo_total=0
nhi_total=0

FOR j=0,nexp-1 DO BEGIN
  img=reform(wdout[*,*,j])
  errimg=reform(errout[*,*,j])
  k=where(errimg EQ missing,nk)
  IF nk NE 0 THEN img[k]=missing
  
  s=size(img,/dim)
  swtch_lo=make_array(s,/byte)
  swtch_hi=make_array(s,/byte)
  swtch_pix=make_array(s,/byte,value=1b)
    
  IF max(img) GT missing THEN BEGIN

    IF nk NE 0 THEN swtch_pix[k]=0b

    h=histogram(img,min=-30,max=50,bin=1,loc=loc)
    getmax=max(smooth(h,5),imax)
    max_var_low=max([loc[imax],6.])

    k=where(img NE missing AND img LE max_var_low*3.)
    s=max([stdev(img[k]),5.])

   ; re-define max_var_low based on s
   ;
    max_var_low=max([loc[imax]+2.0*s,20.])

   ;
   ; This kernel gives the number of non-missing pixels used in the
   ; 3x7 median calculation ("med3").
   ;
    kernel=make_array(3,7,value=1b)
    npix=convol(swtch_pix,kernel,/edge_zero)
   ;
   ; These kernels give the number of pixels above and below the
   ; central pixel in the 1x7 median calculation ("med1").
   ;
    kernel=bytarr(1,7)
    kernel[0,*]=[1b,1b,1b,0b,0b,0b,0b]
    nbelow=convol(swtch_pix,kernel,/edge_zero)
    kernel[0,*]=[0b,0b,0b,0b,1b,1b,1b]
    nabove=convol(swtch_pix,kernel,/edge_zero)

    limit=3.*max_var_low
    med3=fmedian(img,3,7,missing=missing)
    med1=fmedian(img,1,7,missing=missing)

   ;
   ; This is the data structure that gets passed to
   ; EIS_CR_CHECKER. 
   ;
    data={ img: img, $
           med3: med3, $
           med1: med1, $
           npix: npix, $
           nabove: nabove, $
           nbelow: nbelow,$
           limit: limit, $
           max_var_low: max_var_low, $
           sigma: s, $
           missing: missing, $
           swtch_hi: swtch_hi, $
           swtch_lo: swtch_lo }

   ;
   ; This is the first pass on the image. Strong ("hi") and weak
   ; ("lo") CRs are flagged separately in the swtch_lo and swtch_hi
   ; arrays in DATA.
   ;
    eis_cr_checker, data
    swtch_lo=data.swtch_lo
    swtch_hi=data.swtch_hi
    
   ;
   ; The use of kernel and convol here makes the neighbors of the
   ; flagged pixels in swtch_hi to be flagged too.
   ;
    kernel=[[0b,1b,0b],[1b,1b,1b],[0b,1b,0b]]
    touched=convol(swtch_hi,kernel,/edge_zero)
    k=where(touched GT 0,nk)
    IF nk NE 0 THEN swtch_hi[k]=1b

   ;
   ; The following does a second pass on the high intensity CRs to
   ; remove any residuals (which seem to be quite common). The med1
   ; and med3 arrays need to be re-computed to avoid using the CRs
   ; from the first pass.
   ;
    img2=img
    swtch_all=swtch_lo OR swtch_hi
    k=where(swtch_all EQ 1b,nk)
    IF nk NE 0 THEN img2[k] = missing

    data.med3=fmedian(img2,3,7,missing=missing)
    data.med1=fmedian(img2,1,7,missing=missing)
    data.img=img2
    data.swtch_hi=swtch_hi

   ;
   ; This is the second pass on the image, but with the updated MED3
   ; and MED1 arrays
   ;
    eis_cr_checker,data
    swtch_lo=data.swtch_lo
    swtch_hi=data.swtch_hi

   ;
   ; Now combine the lo and hi CR arrays and print statistics.
   ;
    k=where(swtch_lo EQ 1,n1)
    k=where(swtch_hi EQ 1,n2)
    swtch[*,*,j]=swtch_hi OR swtch_lo

   ;
   ; I've commented this out, but I leave it in case of future
   ; bug-checking.
   ;
;    IF keyword_set(verbose) THEN print,'% EIS_PREP_CLEAN_CRS: Exp No. '+trim(j)+', pixels removed: '+trim(n1)+', '+trim(n2)+'   max_var_low: '+trim(string(max_var_low,format='(f8.1)'))+'  stdev: '+trim(string(s,format='(f8.1)'))

    nlo_total=nlo_total+n1
    nhi_total=nhi_total+n2
  ENDIF 
ENDFOR

IF keyword_set(verbose) THEN BEGIN
  print,'% EIS_PREP_CLEAN_CRS: number of low intensity CRs = '+trim(nlo_total)
  print,'% EIS_PREP_CLEAN_CRS: number of high intensity CRs = '+trim(nhi_total)
ENDIF 

k=where(swtch EQ 1b,nk)
IF nk NE 0 THEN errout[k]=missing

;
; Now fill the warm pixels with interpolated values, or set them to
; missing if /nofill set.
;
IF keyword_set(nofill) THEN BEGIN
ENDIF ELSE BEGIN
  FOR i=0,nexp-1 DO BEGIN
    swtch_img=reform(swtch[*,*,i])
    index=where(swtch_img EQ 1,ni)
    IF ni NE 0 THEN BEGIN
      img=reform(wdout[*,*,i])
      wdout[*,*,i]=eis_prep_med_fill(img,index=index,missing=missing)
    ENDIF 
  ENDFOR 
ENDELSE 


IF keyword_set(swap_xy) THEN BEGIN
  swtch=transpose(temporary(swtch),[0,2,1])
  wdout=transpose(temporary(wdout),[0,2,1])
  errout=transpose(temporary(errout),[0,2,1])
ENDIF 

err=temporary(errout)
wd=temporary(wdout)

END