PRO xrt_clean_ro, image_in, index_in, dark, image_out, grade, model,   $
                  hist=hist, clean_type=clean_type, $
                  bsub_type=bsub_type, dark_type=dark_type, nsigma=nsigma,  $
                  nmed=nmed, stop_bleed=stop_bleed, verbose=verbose, $ 
                  quiet=quiet, qabort=qabort, qstop=qstop

; =========================================================================
;+
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;       XRT_CLEAN_RO
;
; CATEGORY:
;       Data calibration
;
; PURPOSE:
;       Remove periodic and other "background" signals introduced into 
;       the data by the read-out electronics.
;
; CALLING SEQUENCE:
;       XRT_CLEAN_RO, image_in, index_in, dark, image_out  [,grade]  
;                    [,model] ,hist=hist [,clean_type=clean_type] 
;                    [,bsub_type=bsub_type] [,dark_type=dark_type] 
;                    [,nsigma=nsigma] [,nmed=nmed] [,/stop_bleed] 
;                    [,/verbose] [,/quiet] [,qabort=qabort] [,/qstop]
;
; INPUTS:
;       IMAGE_IN  - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image containing read-out signals.
;       INDEX_IN  - [Mandatory] index structure associated with IMAGE_IN
;          DARK   - [Mandatory] (medianed dark) to be used for either
;                   dark subtraction or background level adjustment. 
;                   Typically generated by <xrt_med_dark.pro>.
;                   (may be non-existant, in which case it is generated)
;          GRADE  - [Optional] a constant indicating the saturation level
;                     (default = 2500 DN)
;                    OR 
;                   a 2-dim byte array, [Nx,Ny] containing 
;                    (at a minimum)
;                    0 if pixel is OK
;                    1 if pixel is saturated (>2500 DN)
;                    (note: overwritten on output)
;
; KEYWORDS:
;           HIST  - [Mandatory] (string array) Contains information
;                   on type of cleaning done and keyword values
;                   for later use in updating header information
;      CLEAN_TYPE - [Optional] Type of Fourier cleaning to use:
;                   0 = prefilter on semi-fixed streaks, then remove   
;                       remaining Fourier streaks, stars, & smudges 
;                       [default] 
;                   1 = remove semi-fixed streaks only
;                   2 = remove Fourier streaks, stars, & smudges
;                   3 = don't Fourier clean!
;       BSUB_TYPE - [Optional] Type of background subtraction to use:
;                   0 = remove Nyquist ringing and large-scale background
;                       "ramp" (dark) [default]  
;                   1 = remove Nyquist ringing only
;                   2 = remove large-scale background "ramp" (dark) only 
;                   3 = no (model) background subraction!
;       DARK_TYPE - [Optional] Type of dark subtraction to use:
;                   0 = Use model dark with median adjusted to median of 
;                       cleaned darks nearby in time (if possible) [default]
;                   1 = Use median of cleaned darks nearby in time 
;                       (if possible)  
;                   2 = Use model darks without median level adjustment
;                   (Note: if dark_type=1, bsub_type = 3 is used)
;         NSIGMA  - [Optional] Number of standard deviations
;                   beyond which an FFT pixel is
;                   considered bad in "Fourier star" removal. 
;                   Default is 4.5
;                   (Note: much below 4 not recommended)
;         NMED    - [Optional] Number of standard deviations
;                   above smoothed central minimum background
;                   to begin shielding (presumed) data from correction
;                   Default is 3.5
;                   (Note: outside range of 2.0 to 4.5 not recommended)
;     /STOP_BLEED - [Optional] (Boolean) Set = 1 if pixels corrected for
;                   saturation "bloom/bleed" at the edge of saturated
;                   region(s) are to be retained (otherwise they are
;                   reverted to their input values)
;       /QUIET    - [Optional] (Boolean) Set for no messages or
;                   errors. Default is some messages and all errors.
;       /VERBOSE  - [Optional] (Boolean) Set for all messages, errors, and
;                   intermediate data listings. Suppresses quiet.
;       /QSTOP    - [Optional] (Boolean) For debugging.
;
; OUTPUTS:
;       IMAGE_OUT - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image with the model of read-out signals removed.
;          GRADE  - [Optional] a 2-dim byte array (1x1 binning case) or
;                    a 2-dim double array (other binnings), [Nx,Ny] containing 
;                    the sum of
;                    0 if pixel is OK
;                    1 if pixel is saturated (default >2500 DN)
;                    2 if pixel is a victim of saturation "bloom/bleed"
;                    4 if pixel is in a contamination spot
;                    8 if pixel is in a dust speck
;                   16 if pixel is "hot" 
;                   PLUS (optionally, in case of binning, if the /grade_frac
;                   flag is set) an encoding (in fractional form).
;       MODEL     - [Optional] (2-dim float array, [Nx,Ny])
;                   The model of read-out signals which has been removed
;                   (subtracted).
;       QABORT    - [Optional] (Boolean) Indicates that the program exited
;                   gracefully without completing. (Might be useful
;                   for calling programs.)
;                   0: Program ran to completion.
;                   1: Program aborted before completion.
;
; EXAMPLES:
;
;       Basic usage: 
; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, hist=hist
;
;      You want the the pixel grade and read-out/dark model: 
; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, grade, model, $
;         hist=hist
;
;       You want to correct saturation bleed-affected pixels, and 
;          pixel grades: 
; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, grade, $ 
;          hist=hist, stop_bleed=1 
;
;       You want pixel grades, but no Fourier cleaning:
; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, grade, $ 
;          hist=hist, clean_type=3 
;
;       You want more aggressive Fourier cleaning, with median dark
;         subtraction:
; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, $ 
;          hist=hist, nsigma=4, nmed=4, dark_type=1
;
;
;
; COMMON BLOCKS:
;       None.
;
; NOTES:
;       1) An XRT image shows a prevelent sinusoidal ringing in x (width) 
;        at the nyquist frequency, with a peak-to-peak amplitude of 
;        ~2.6 DN due to different bias levels set in odd vs. even columns.  
;        There is also a low-level "ski-ramp" in the y (height)
;        direction with a shape which can be approximated by an exponential
;        decrease (amplitude ~ 4.3 DN, e-folding width ~185 pixels) + a 
;        weak linear increase on a base of ~42 DN (for 1x1 binning).  
;        (This is best seen in dark frames).  The Fourier power spectrum of 
;        an XRT dark shows many 
;        odd features, including "streaks" (narrow ridges of power spanning
;        all y, many are semi-fixed in location), "stars"  (2-D Voigt 
;        profiles) and "smudges" (truncated streaks).  The latter two 
;        features are typically variable in y position; all vary
;        in amplitude. These features are present (though less clearly 
;        visible) in the Fourier transform of data as well.  This 
;        procedure attempts to correct for these read-out signals.
;       2) The Fourier filtering can take a while on large images
;        (~15s on a 2048x2048 using a MacPro with dual 3GHz Xeon chips). 
;        Its effects can be dramatic or minimal, but are typically 
;        confined to low count levels.  
;        Thus, for certain types of images and certain applications, the 
;        user may want to skip this part of the cleaning (set clean_type=3).  
;        See xrt_fourier_vacuuum.pro notes for more details
;        on how to use and modify the Fourier cleaning.
;       3) Some Fourier patterns may not be fully removed in images with one
;        or more smaller dimensions (<=128).  Fourier patterns are not
;        removed if a large part of the image (>45-80% depending on
;        image size) is saturated.
;
; MODIFICATION HISTORY:
;
progver = 'v2007.Feb.14' ; --- (SS, MW) Written by SS. Cleaned up by MW.
progver = 'v2007.Mar.06' ; --- (SS, MW) Modified by SS to include 
                         ; --- large-scale background removal, reorganize 
                         ; --- using function calls . Cleaned up by MW.
progver = 'v2007.Apr.23' ; --- (SS, MW) Modified by SS to include a first 
                         ; --- stab at Fourier read-out signal removal. 
                         ; --- Cleaned up by MW.
progver = 'v2007.May.18' ; --- (SS) Modified by SS to include
                         ; --- "patch-over" of saturated regions to
                         ; --- improve Fourier read-out signal 
                         ; --- removal, plus a few other minor fixes.
progver = 'v2007.Nov.27' ; --- (SS) Added comments/history update about 
                         ; --- possible poorer performance for small images.
progver = 'v2007.Nov.28' ; --- (SS) Added fix for cases where images are
                         ; --- significantly saturated.
progver = 'v2008.May.18' ; --- (SSaar) Added an improved large-scale 
                        ; --- background correction (see <lsback_away.pro>)
                        ; --- including corrections for binning, CCD 
                        ; --- temperature, and exposure time.
progver = 'v2009.Jun.12' ; --- (SSaar) Added an improved dark subtraction
                        ; --- adjusting model results with median of darks 
                        ; --- observed nearby in time.  Added generation
                        ; --- of initial pixel "grade" map.  Added options for 
                        ; --- different kinds of dark subtraction. Made
                        ; --- some minor corrections.
progver = 'v2009.Jul.08' ; --- (SSaar) Fixed a bug in dark zero level reset
progver = 'v2009.Nov.10' ; --- (SSaar) Fixed a bug when i_bleed does not exist
;                        ; ---- or = -1 
progver = 'v2010.Dec.02' ; --- (SSaar) Fixed a bug when bin=4 and improved
;                        ; --- correction for dark_type=1
progver = 'v2012.Feb.15' ; --- (SSaar) Fixed a bug when i_bleed does 
                         ;  not exist and clean_type=3
progver = 'v2013.Feb.15' ;--- (SSaar (SAO)) Fixed problem arising when
;                        ;---  bad images are mislabeled as darks.
;                        ;---  Simplified the code structure somewhat; 
;                        ;---  improved reuse of already processed median darks
;                        ;--- (should speed code somewhat).
progver = 'v2014.Feb.06' ;--- (JSlavin (SAO)) Pass along quiet flag to
                         ;--- xrt_med_dark
;
;-
; =========================================================================
;=== Check inputs & initial setup =====================================

;===== Set some keyword Booleans, plus miscellaneous.
  q_vb = keyword_set(verbose)
  q_qt = keyword_set(quiet) and (not q_vb)
  qstop        = keyword_set(qstop)
  qabort       = 0B
  if not keyword_set(clean_type) then clean_type=0  ; default to full Fourier
  if not keyword_set(bsub_type) then bsub_type=0  ; default to full background
                                                ; subtraction
  if not keyword_set(dark_type) then dark_type=0  ; default to median adjusted 
                                                  ; dark model 
  if qstop eq 1 then debug=1 else debug=0

  if not keyword_set(nsigma) then nsigma=4.5 ; threshold sigma to flag spike
  if not keyword_set(nmed) then nmed=3.5 ; shield data above nmed*central_median
  if not keyword_set(stop_bleed) then stop_bleed=0 ; do not replace bleed pixels

  if n_elements(dark) lt 1 then $             ; calc medianed dark
    dark = xrt_med_dark(index_in, ddates=dates, darkff=darkff, nday=15, $
                           quiet=q_qt) 
;stop
;===== Thresholds for "small" image size, maximum saturated pixel % 
  sat_level = 2500.
  if n_elements(grade) eq 0 then grade=sat_level
  minsize = 128.
  maxsat = 0.80

;===== Check inputs.

  if q_vb then print, 'XRT_CLEAN_RO: Checking input parameters...'

  ;=== Check number of parameters.
  if (n_params() ge 3) then begin
      ; ===== OK
  endif else begin

    if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires >=3 ' $
                                    + 'parameters. Aborting...'
    qabort = 1B
    return
  endelse 

  ;=== Check size and type of input array.
  dim_check = (size(image_in,/n_dim) eq 2)
  typ_check = (datatype(image_in) eq 'FLO') or (datatype(image_in) eq 'DOU')
  if (dim_check and typ_check) then begin
    ;=== OK
  endif else begin
    if not dim_check then begin
      if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires a 2D ' $
                                 + 'image array as input. Aborting...'
      qabort = 1B
    endif
    if not typ_check then begin
      if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires that ' $
                           + 'the image is floating type. Aborting...'
      qabort = 1B
    endif
    return
  endelse

;===== check clean_type
  if (clean_type lt 0) or (clean_type gt 3) or $ 
                           (clean_type ne fix(clean_type)) then begin
       if (not q_qt) then box_message, 'XRT_CLEAN_RO:  Requires that ' $
         + 'clean_type = 0, 1, 2, or 3. Defaulting to clean_type = 0'
      clean_type=0
  endif

;===== check bsub_type
  if (bsub_type lt 0) or (bsub_type gt 3) or $ 
                           (bsub_type ne fix(bsub_type)) then begin
       if (not q_qt) then box_message, 'XRT_CLEAN_RO:  Requires that ' $
         + 'bsub_type = 0, 1, 2, or 3. Defaulting to bsub_type = 0'
      bsub_type=0
  endif

  if q_vb then print, '                                       ...OK'

;===== Get dimensions of image.
  s = size(image_in)
  n1 = s(1)
  n2 = s(2)
  x1 = findgen(n1)
  bin = index_in.chip_sum

;===== If Fourier cleaning is used, temporarily smooth over saturated
;=====   areas (as needed) and keep record of altered "bloom/bleed" pixels 
;=====   If number of saturated pixels > 80% (128x128) to 45% (2048x2048)
;=====   of the total, default to no Fourier cleaning (clean_type=3)
;=====   and add note to History later

  satflag=0
  i_bleed=-1         ; covers case of sat exist but sat_fudge not used 
  if clean_type ne 3 then begin
     if n_elements(grade) gt 1 then nsat = total(grade eq 1) else $ 
        nsat = total(image_in ge grade)
     if nsat gt (n1*n2)*maxsat/(n1*n2/minsize^2)^0.1 then satflag=1
     if satflag eq 0 then begin
           if q_vb then print,'Fudging saturated pixels (if any)...' 
           saturation_fudge,image_in,grade,image,i_bleed 
     endif else begin
           if q_vb then print,'Defaulting to no Fourier cleaning - too' $
                                 + ' many saturated pixels!'
           clean_type=3
           image = image_in
     endelse
  endif else image = image_in

  if q_vb then begin
      case 1 of
         (bsub_type eq 1): txtv='  Remove Nyquist ringing only'
         (bsub_type eq 2): txtv='  Remove large-scale ramp only'
         (bsub_type eq 3): txtv='  No Nyquist or ramp subtraction'
         else: txtv='  [Default] Remove Nyquist ringing and large-scale ramp'
      end
      print, 'XRT_CLEAN_RO:  beginning Nyquist and ramp subtraction; '
      print,txtv
  endif
 
;print,'ct: ',clean_type,'bt: ',bsub_type,'dt: ',dark_type 

  dflag=(n_elements(dark) gt 1)       ; 0 if no valid dark
  if dflag eq 0 then mdark = 0. else mdark = median(dark)
  txt='XRT_CLEAN_RO: '
  if bsub_type ne 3 then begin      ; (if so, then dark_type <> 1)
     if bsub_type ne 2 then begin     ; if not "ramp only"
;======  Remove Nyquist Frequency ringing
        image_out1 = no_nyquist(image, index_in, fit1)
        txt=txt+'  Remove Nyquist ringing'
     endif else image_out1 = image
     if bsub_type ne 1 then begin          ; if not "Nyquist only"
;===== Remove Large-Scale component of background 
;=====        (Pedestal+Read-out+Dark) 
        txt=txt+'  Remove large-scale background ramp '
        if dark_type eq 0 then begin  ; default "optimum" case: subtract model 
                                ; with median adjusted to medianed dark, if 
                                ; at all possible (else just median or model)
           if (bin ne 4 or (bin eq 4 and dflag eq 0)) then begin
              image_out2 = lsback_away(image_out1, index_in, fit2)
              image_out2 = image_out2 -mdark +median(fit2)*dflag 
              fit2 = fit2 +mdark -median(fit2)*dflag 
              if dflag eq 0 then txt=txt+' via model dark.' else $
                 txt=txt+' via median adjusted model dark [default].'
           endif else begin     ;  use median darks where possible if bin=4
              image_out2 = image_out1 - dark
              fit2 = dark
              txt=txt+' via medianed dark.'
           end
        endif else begin             ; dark_type=3, subtract model dark
           image_out2 = lsback_away(image_out1, index_in, fit2)
           txt=txt+' via model dark.'
        end
     endif 
  endif else begin 
     if dark_type eq 1 then begin   ; subtract median dark if at all possible
        image_out1 = no_nyquist(image, index_in, fit1)
        if total(mdark) ne 0 then begin
           image_out2 = image_out1 - dark
           fit2 = dark
           txt=txt+' via medianed dark.'
        endif else begin        ; if no good darks, use model
           image_out2 = lsback_away(image_out1, index_in, fit2)
           txt=txt+' via model dark [not enough good darks].'
        end
     endif else begin
        txt=txt + '  NO Nyquist or large-scale ramp removal.
        image_out2 = image
     end
  end 
  if bsub_type eq 0 then txt=txt+'  [Default]'
;===== construct history update
  hist=strarr(3)
  hist(0)=txt


  if q_vb then begin
      print, 'XRT_CLEAN_RO:  beginning Fourier cleaning '
      case 1 of
         (clean_type eq 1): txt='  Remove semi-fixed streaks only'
         (clean_type eq 2): txt='  Remove Fourier streaks, stars, & smudges'
         (clean_type eq 3): txt='  No Fourier cleaning.'
         else: txt='  [Default] Prefilter on semi-fixed streaks, then ' $ 
                    + ' remove remaining Fourier streaks, stars, & smudges'
      end
      if satflag eq 1 then txt = txt + '  Too many saturated pixels! '
      print, 'clean_type = ',clean_type,' '+txt
      if nsigma eq 4.5 then txt=' [default] ' else txt=''
      print,' nsigma = ',nsigma,txt
      if nmed eq 3.5 then txt=' [default] ' else txt=''
      print,' nmed = ',nmed,txt
  endif	


;=====  Vacuum up naughty read-out features in Fourier space
  if clean_type ne 3 then begin  
      xrt_fourier_vacuum,image_out2,image_out, $
         clean_type=clean_type,nsig=nsig,nmed=nmed,debug=debug

;===== if stop_bleed=0, change "bleed/bloom" boundary pixels to 
;=====        original values with nyquist and ramp correction only
      if ((stop_bleed eq 0) and (i_bleed[0] ne -1)) then begin 
         case 1 of
            (bsub_type eq 0): corr = fit1(i_bleed)+fit2(i_bleed)
            (bsub_type eq 1): corr = fit1(i_bleed)
            (bsub_type eq 2): corr = fit2(i_bleed)
            else: corr = 0.
         end
         image_out(i_bleed)=image_in(i_bleed) - corr
      endif

  endif else image_out = image_out2

;===== add to history update
  case 1 of
      (clean_type eq 1): txt=' Semi-fixed Fourier streak preclean'
      (clean_type eq 2): txt=' Fourier cleaning (no streak preclean)'
      (clean_type eq 3): txt=' No Fourier cleaning. '
       else:             txt=' Full Fourier preclean/clean [default]'
  end
  if satflag eq 1 then txt = txt + '  Too many saturated pixels! '
  hist(1)='XRT_CLEAN_RO: clean_type = ' $ 
    + strtrim(string(fix(clean_type)),2) + txt
  if clean_type ne 3 then begin 
     hist(2)= 'XRT_CLEAN_RO: nmed = ' + strtrim(string(nmed),2)
     if nmed eq 3.5 then hist(2) = hist(2) + ' [default]'
     if clean_type ne 1 then begin
        hist(2)= hist(2) + '  nsigma = ' + strtrim(string(nsigma),2)
        if nsigma eq 4.5 then hist(2) = hist(2) + ' [default]'
     endif
  endif else hist=hist(0:1)
  if ((n1 le minsize) or (n2 le minsize)) and ((clean_type eq 0) or $
    (clean_type eq 2)) then hist = [hist,'XRT_CLEAN_RO: small image '+ $
     'dimension - some Fourier noise patterns may remain uncorrected']
  if stop_bleed eq 1 and satflag eq 0 then hist = $
     [hist,'XRT_CLEAN_RO: saturation-affected boundary pixels corrected ']
  hist=['XRT_CLEAN_RO: Running '+progver,hist]  

;=====  reset saturation level to input value
if n_elements(grade) eq 1 then isat=where(image_in ge grade,nsat) $ 
    else isat=where(grade eq 1,nsat)
if nsat gt 0 then image_out(isat) = image_in(isat)

;===== generate fit array 
  model = image_in - image_out

;====  generate initial grade array (saturation and bleed only)
  grade=bytarr(n1,n2)
  if nsat gt 0 then begin
      grade(isat) = 1B   
      if n_elements(i_bleed) gt 0 and i_bleed(0) ne -1 then grade(i_bleed) = 2B
                                    ; checks i_bleed is there and not = -1 
  endif

END ;======================================================================