FUNCTION eis_prep_med_fill, img, err, index=index, quiet=quiet, missing=missing

;+
; NAME:
;     EIS_PREP_MED_FILL
;
; PURPOSE:
;     This routine contains the code for filling in missing data using
;     a median of neighbors. The procedure is taken directly from
;     eis_prep. 
;
; CATEGORY:
;     Hinode/EIS; calibration.
;
; CALLING SEQUENCE:
;     Result = EIS_PREP_MED_FILL( Img )
;
; INPUTS:
;     Img:   A 2D image assumed to be an EIS intensity array. If
;            neither ERR or INDEX are input then missing data
;            should be marked with the value MISSING. If INDEX is
;            input, then IMG[INDEX] will be set to MISSING. If ERR is
;            input (and INDEX isn't) then INDEX is set to where
;            the pixels where ERR is set to MISSING.
;
; OPTIONAL INPUTS:
;     Index: A 1D index array giving the locations of the pixels to be
;            filled. 
;     Err:   The 2D error array corresponding to IMG.
;     Missing:  The value for missing data. If not set, then the value
;               -100 is assumed.
;	
; KEYWORD PARAMETERS:
;     QUIET:  If set, then no information messages will be printed to screen.
;
; OUTPUTS:
;     A 2D image of the same size as IMG, but with either missing
;     pixels (if ERR specified) or selected pixels (if INDEX
;     specified) replaced be neighboring pixel values.
;
;     If there's a problem then IMG is returned to the user.
;
; EXAMPLE:
;     Fill all missing pixels:
;      IDL> imgout=eis_prep_med_fill(img,err)
;
;     Fill only pixels stored in index:
;      IDL> imgout=eis_prep_med_fill(img,index=index)
;
; PROGRAMMING NOTES:
;     When used with EIS_PREP you should use INDEX to flag missing
;     pixels, rather than ERR (or INT). This is because EIS_PREP
;     progressively fills the different types of missing data (saturated,
;     warm pixels, cosmic rays, etc.) in INT, whereas ERR accumulates all
;     missing data. INDEX can be used to flag each missing data type
;     as they're processed so that we don't try to repeatedly fill all
;     of the accumulated missing pixels in ERR.
;
; MODIFICATION HISTORY:
;     Ver.1, 20-Sep-2019, Peter Young
;        Code adapted from eis_prep.
;-

IF n_params() LT 1 THEN BEGIN
  print,'Use:  IDL> imgout = eis_prep_med_fill( img [, err, index=index, missing=missing, /quiet ])'
  print,'   either ERR or INDEX must be input'
  return,img
ENDIF 

IF n_elements(err) EQ 0 AND n_elements(index) EQ 0 THEN BEGIN
  print,'% EIS_PREP_MED_FILL: please specify at least one of ERR or INDEX. Returning...'
  return,img
ENDIF

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

s=size(img)
IF s[0] NE 2 THEN BEGIN
  print,'% EIS_PREP_MED_FILL: the input IMG must be a 2D array. Returning...'
  return,img
ENDIF 

imgout=img

IF max(img) EQ missing THEN BEGIN
  IF NOT keyword_set(quiet) THEN print,'% EIS_PREP_MED_FILL: the image contains only missing data. Returning...'
  return,img
ENDIF 


IF n_elements(index) EQ 0 THEN index=where(err EQ missing,ni) ELSE ni=n_elements(index)
IF ni NE 0 THEN imgout[index]=missing ELSE return,imgout
  
  
ct=ni
xbox=1 & ybox=3
;
; The procedure is to fill missing data with the median of the two
; neighbors in the Y-direction. If the two neighbors are also missing
; then the routine just repeats until good data is moved into the
; neighboring pixels. If an entire column is missing then
; interpolation is not possible and the procedure will stall at
; ct=ct_old. 
;
REPEAT BEGIN
  ct_old=ct
  newmed = fmedian(imgout,xbox,ybox,missing=missing,/only_missing)
  imgout[index] = newmed[index]
  newmiss = where(newmed[index] eq missing,ct)
END UNTIL (ct eq 0) or (ct ge ct_old)

return,imgout

END