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

;+
; NAME:
;     EIS_RESIDUAL_WPS
;
; PURPOSE:
;     Identify residual warm pixels in 3D EIS data arrays. Note this
;     is intended to be called from within EIS_PREP and
;     generally should not be called separately. An EIS Software Note
;     gives more details about the method.
;
; CATEGORY:
;     Hinode; EIS; calibration; warm pixels.
;
; CALLING SEQUENCE:
;     EIS_PREP_RESIDUAL_WPS, 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
;
; PROGRAMMING NOTES:
;     This routine should be called from eis_prep after the hot 
;     and warm pixels have been flagged, and prior to cosmic rays
;     being flagged.
;
;     To perform tests on the routine it is recommended that you run
;     eis_prep with the /noabs and /nocr keywords to create a level-1
;     file. Then load a data window using eis_getwindata and run
;     eis_clean_windata as follows:
;
;     IDL> wd=eis_getwindata(l1name,wvl)
;     IDL> wdout=eis_clean_windata(wd,swtch=swtch,/nocr)
;
;     The pixels that have been flagged can be browsed with:
;
;     IDL> eis_spike_viewer,wd,wdout,swtch
;
; EXAMPLE:
;     See the routines EIS_PREP and EIS_CLEAN_WINDATA for how to use
;     this routine. 
;
; MODIFICATION HISTORY:
;     Ver.1, 5-Nov-2019, Peter Young
;     Ver.2, 17-Jun-2020, Peter Young
;       Fixed bug when setting flags for very strong pixels. 
;-

IF n_params() LT 2 THEN BEGIN
  print,'Use:  IDL> eis_prep_residual_wps, 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 4 THEN BEGIN
  print,'% EIS_PREP_RESIDUAL_WPS: The input WD has 4 dimensions and is not compatible with this routine. Returning...'
  delvarx,wdout,errout
  return
ENDIF 
IF s[0] EQ 2 THEN nexp=1 ELSE nexp=s[3]
s=size(wdout,/dim)
swtch_wp=make_array(s,/byte)
swtch=swtch_wp
swtch_miss=wdout EQ missing

;
; 'Free' parameters (they are not free as I've hard-coded them!
; but if you want to tune the method, these are the numbers to use.
;
min_exp=10      ; minimum no. of exposures the method needs
min_dn=6        ; minimum DN level for identifying weak WPs
strong_scl=1.5  ; scale parameter for strong pixel check.
exp_frac=0.50   ; fraction of exposures required to flag WP


;
IF nexp LT min_exp THEN BEGIN
  IF keyword_set(verbose) THEN print,'% EIS_PREP_RESIDUAL_WPS: Data window has less than 10 exposures. Residual WPs can not be flagged. Returning...'
  delvarx,wdout,errout,swtch_wp,swtch_miss
  return
ENDIF 


;
; For identifying potential warm pixels in a 2D image, I do:
;   1. Compare with the median of a 3x7 pixel block. The pixel must
;      satisfy:
;      a. pix-med > max_var_low
;      b. pix < limit
;      c. the 3x7 block must have 7 or more non-missing pixels.
;   2. Compare with the median of a 1x5 pixel block. The pixel must
;      satisfy:
;      a. pix-med > simg*1.5
;      b. pix > limit
;      c. the 1x5 block must have 5 or more non-missing pixels.
;      d. pix must be < 1000 DN
;   3. Very bright pixels should satsify:
;      a. pix-med > simg*10
;      b. img > 1000
;      c. the 1x5 block (see above) must have 5 or more non-missing
;         pixels.
;
; If a pixel satisfies any of the above criteria, then I require that
; the same pixel is flagged in 60% of the exposures or more than 10
; pixels (whichever is the larger number of pixels) to be flagged as a
; residual warm pixel.
;
; The parameters MAX_VAR_LOW, LIMIT and SIMG are computed for each
; individual exposure. They are not fixed parameters.
;
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_img=make_array(s,/byte)
  swtch_pix=make_array(s,/byte,value=1b)
    
  IF max(img) EQ missing THEN CONTINUE

 ;
 ; 'npix' is the number of pixels that enter into the median
 ; calculation, and is determined from swtch_pix which is an array
 ; of 1's, except for the positions of missing pixels, which are
 ; 0's. Using the convol function on this array allows npix
 ; to be determined.
 ;
  swtch_pix[k]=0b
  kernel=make_array(3,7,value=1b)
  npix=convol(swtch_pix,kernel,/edge_zero)
    
 ;
 ; When doing the median check on 1x7 blocks (see below) we want to
 ; make sure there are pixels on both sides of the pixel that are
 ; being used for the median calculation. For this I compute
 ; 'nbelow' and 'nabove' using suitable kernels and the convol
 ; function.
 ;
  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)  

 ;
 ; I take a histogram of weak count pixels (-30 to 50 DN). The peak of
 ; this histogram is taken as the noise level (loc[imax]). Thus if a
 ; pixel is loc[imax] above the median, then there's a good
 ; chance it could be a warm pixel. loc[imax] can't be too low,
 ; however, so I force it to be > 6.
 ;
  h=histogram(img,min=-30,max=50,bin=1,loc=loc)
  getmax=max(smooth(h,5),imax)
 ;
 ; Define parameters.
 ;
  max_var_low=max([loc[imax],min_dn])
  limit=(max_var_low/1.5)^2
  simg=sqrt(img>0)>max_var_low

 ;
 ; This is check on weak pixels [ < limit DN ]
 ; 
  med3=fmedian(img,3,7,missing=missing)
  k=where( (img-med3) GT max_var_low AND (med3 NE missing) AND (img LE limit) AND (npix GE 7), nk)
  IF nk NE 0 THEN swtch_img[k]=1b

 ;
 ; This is check on strong pixels [ > limit DN ]
 ;
 ; Here I'm checking the 1x7 median (not 3x7).
 ; Initially I required the pixel to be sqrt(N) larger than the median
 ; since is the noise level. I scaled this up empirically by looking
 ; at the flagged pixels.
 ; 
  med1=fmedian(img,1,7,missing=missing)
  k=where( ((img-med1) GT simg*strong_scl) AND (img GT limit) AND (med1 NE missing) AND (nbelow GE 1) AND (nabove GE 1) AND (img LE 1000), nk)
  IF nk NE 0 THEN swtch_img[k]=1b

 ;
 ; This is check on very strong pixels [ > 1000 DN ]
 ; 
 ; I require the pixel
 ; to be >10sigma above the median, so this will be a very rare
 ; event (an example is 20-jan-2007 22:32, 195 window, pix
 ; [40,171]. Most of these very strong pixels should have been removed
 ; by the hot or warm pixel maps.
 ;
  k=where( ((img-med3) GT simg*10) AND (img GT limit) AND (med3 NE missing) AND (npix GE 3) AND (img GT 1000), nk)
  IF nk NE 0 THEN swtch_img[k]=1b
    
  swtch_wp[*,*,j]=swtch_img

ENDFOR

;
; 'swtch_sumexp' gives the no. of flagged pixels summed over the
;   exposure direction.
; 'swtch_summiss' gives the no. of missing pixels summed over the
;   exposure direction.
;
; I require:
;   - must be less than 50% missing pixels in the exposure direction
;   - no. of non-missing pixels must be > min_exp
;   - the fraction of flagged vs. non-missing pixels must be > 60%.
;
; Exposure pixels that satisfy these criteria are flagged as residual
; warm pixels.
;
swtch_sumexp=total(swtch_wp,3)
swtch_summiss=total(swtch_miss,3)
i=where( (nexp-swtch_summiss) GE 0.5*float(nexp) AND (nexp-swtch_summiss) GE min_exp AND (float(swtch_sumexp)/float(nexp-swtch_summiss) GE exp_frac),ni)

perc=float(ni)/float(n_elements(swtch_sumexp))*100
IF keyword_set(verbose) THEN print,'% EIS_PREP_RESIDUAL_WPS:  No. of residual warm pixels = '+trim(ni)+' ['+trim(string(perc,format='(f6.1,"%")'))+']'

;
; Fill the SWTCH output array and update ERR with the new missing values.
;
IF ni GT 0 THEN BEGIN 
  FOR j=0,nexp-1 DO BEGIN
    swtch_img=reform(swtch[*,*,j])
    swtch_img[i]=1
    swtch[*,*,j]=swtch_img
  ENDFOR
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
  wdout[k]=missing
ENDIF ELSE BEGIN
  IF ni GT 0 THEN BEGIN
    index=where(swtch_img EQ 1)
    FOR i=0,nexp-1 DO wdout[*,*,i]=eis_prep_med_fill(reform(wdout[*,*,i]),index=index,missing=missing)
  ENDIF 
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