PRO xrt_clean_ro, image_in, index_in, sat, image_out, i_bleed, model,   $
                  hist=hist, clean_type=clean_type, bsub_type=bsub_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, sat, image_out [,i_bleed]
;                     [,model] [,hist=hist] [,clean_type=clean_type] 
;                    [,bsub_type=bsub_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
;            SAT  - [Mandatory] can be either
;                   a 2-dim byte array, [Nx,Ny] containing 1 if pixel
;                   the pixel is saturated, or 0 otherwise
;                   or
;                   scalar value (floating) of the saturation level
;
; 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" [default] 
;                   1 = remove Nyquist ringing only
;                   2 = remove large-scale background "ramp" only
;                   3 = no background sutraction!
;         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.
;       I_BLEED   - [Optional] (1-dim long array, [Nbleed])
;                   Index of pixels where saturation "bloom/bleed"
;                   has been corrected (if stop_bleed=1) or could a
;                   problem, but has been left uncorrected. 
;       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, sat, image_out, hist=hist
;
;       You want the indicies of saturation bleed-affected pixels: 
; IDL> xrt_clean_ro, image_in, index_in, sat, image_out, i_bleed, hist=hist
;
;       You want to correct saturation bleed-affected pixels and 
;       record their indicies: 
; IDL> xrt_clean_ro, image_in, index_in, sat, image_out, i_bleed, $ 
;        hist=hist, stop_bleed=1
;
;       You want the background model for some reason:
; IDL> xrt_clean_ro, image_in, index_in, sat, image_out, i_bleed, model, $ 
;         hist=hist
;
; 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.  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.  (This is best seen 
;        in dark frames).  The Fourier transform 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.
;
;
;-
; =========================================================================
;=== 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 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

;===== Thresholds for "small" image size, maximum saturated pixel % 
  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)

;===== 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
  if clean_type ne 3 then begin
     if n_elements(sat) gt 1 then nsat = total(sat eq 1) else $ 
        nsat = total(image_in ge sat)
     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,sat,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): txt='  Remove Nyquist ringing only'
         (bsub_type eq 2): txt='  Remove large-scale ramp only'
         (bsub_type eq 3): txt='  No Nyquist or ramp subtraction'
         else: txt='  [Default] Remove Nyquist ringing and large-scale ramp'
      end
      print, 'XRT_CLEAN_RO:  beginning Nyquist and ramp subtraction; '
      print,txt
  endif
  

  txt='XRT_CLEAN_RO: '
  if bsub_type ne 3 then begin
     if bsub_type ne 2 then begin
;======  Remove Nyquist Frequency ringing
        image_out1 = no_nyquist(image, fit1)
        txt=txt+'  Remove Nyquist ringing'
     endif
     if bsub_type ne 1 then begin
;===== Remove Large-Scale component of background 
;=====        (Pedestal+Read-out+Dark) 
        image_out2 = lsback_away(image_out1, index_in, fit2)
        txt=txt+'  Remove large-scale background ramp '
     endif
  endif else txt=txt + '  NO Nyquist or large-scale ramp removal.
  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 $ 
        image_out(i_bleed)=image_in(i_bleed)-fit1(i_bleed)-fit2(i_bleed) 

  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(sat) eq 1 then isat=where(image_in ge sat) $ 
    else isat=where(sat eq 1)
if isat(0) ne -1 then image_out(isat) = image_in(isat)

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


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