; Name:
; refsub_h4rg.pro
;
; PURPOSE:
; This program performs a reference pixel correction
; on Rockwell HAWAII-4RG detector data when read out using 32
; outputs.
;
; CALLING SEQUENCE:
;       refsub_h4rg, input_data, output_data, scaName
;
; INPUTS:
;       input_data      - A 4096x4096 array containing the image data.
;       output_data     - A 4096x4096 fltarr for the result
;     scaName       - The name of the device that generated the input data
;
; KEYWORD PARAMETERS:
;       fwindow       - The size of the filter window.
;       order         - Filter order, see the SAVGOL help page
;       degree        - Filter degree, see the SAVGOL help page
;       plot          - Set this =1 to plot the fits
;       simple        - Set this =1 for Maximal Spatial Averaging
;                       (subtract a median reference pixel)
;
;       verbose       - Set this =1 for verbose mode.
;
;
; EXAMPLE
;        IDL> fits_read, 'theImage.fits', myData, myHeader
;        IDL> myNewDate = fltarr(4096, 4096)
;        IDL> refsub_h4rg, myData, myNewData
;
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:
;       Modified:    B.J. Rauscher, IDTL, January 20, 2002
;                    Initial release based on most recent
;                    version of refsub_h1rg
;
;                    B.J. Rauscher, IDTL, February 12, 2003 1/f
;                    noise correlates across all 4 outputs.
;                    Use row reference pixels to make DC
;                    corrections to each output. Use columns to
;                    remove 1/f
;
;             Ernie Morse, IDTL, April 11, 2003 (v2.1)
;             added scaName as input parameter
;             added code block specific to device H2RG-003-4.8mu
;             performed general code cleanup
;
;             Don Figer, RIDL, April 28, 2007
;             modified for H4RG device;
;

; Define a function to return the mean of pixels in reference rows of
; one output. The output number ranges from 0-3.
function get_refrow_mean, data, output

    x_size = 4096./32.; Number of pixels in x-dimension of an output
    _data = fltarr(x_size, 8); Create an array to hold the reference pixel data
    x_start = output * x_size; First column to extract
    x_stop = (output + 1) * x_size - 1; Last column to extract
    _data[0:x_size-1,0:3] = data[x_start:x_stop,0:3]; Get bottom rows of reference pixels
    _data[0:x_size-1,4:7] = data[x_start:x_stop,4092:4095]; Get top rows of refernce pixels

    djs_iterstat, _data, mean = the_mean

    return, the_mean

end


pro refsub_h4rg, input_data, output_data, scaName, fwindow=fwindow, order=order, degree=degree, plot=plot, simple=simple, verbose=verbose, help=help

    ; The reference pixels on the Rockwell HAWAII-4RG SCA are the
    ; outermost 4 pixels all around the chip.

    ; create variables to hold the dimensions of the H2RG SCA (in pixels)
    xsize = 4096
    ysize = 4096

    ; Set defaults
    if (not keyword_set(plot)) then plot = 0 ; Default is no plots
    if (not keyword_set(verbose)) then verbose = 0; Default is quiet mode
    if (not keyword_set(help)) then help = 0; Default does not print help message
    if (not keyword_set(simple)) then simple = 0

    ; Print a help message and exit if requested
    if (help) then $
       stop, 'Usage: refsub_h4rg, input_data, output_data [, fwindow=fwindow, order=order, degree=degree, /plot, /simple, /verbose, /help]'

    ; Adjust all 4 outputs to have the same mean as output #0 in
    ; the reference rows.
    ; Is this needed - We already subtract read 0 which should deal with this? (MWR)
    for output = 1, 3 do begin

       ; Get first and last columns of output to correct
       x_start = output * xsize/4 ; First column of output
       x_stop = (output + 1) * (xsize/4) - 1 ; Last column of output

       ; Offset the output
       input_data[x_start:x_stop,*] = input_data[x_start:x_stop,*] + get_refrow_mean(input_data, output - 1) $
        - get_refrow_mean(input_data, output)

    endfor

    ; We will use a Savitzky-Golay filter to smooth the reference pixel data.
    ; See the SAVGOL procedure's help page for more information on
    ; the following.
    if (not(keyword_set(fwindow))) then fwindow = 64  ; Default uses 64 pixels on
                                                     ; either side of current pixel
                                                     ; for smoothing
    if (not(keyword_set(order))) then order = 0; Use zero to smooth
    if (not(keyword_set(degree))) then degree=3; Three seems reasonable

    ; Extract and average reference pixels
    if (verbose) then $
       print, 'Extracting averaged reference columns on left and right hand sides of SCA'
    row = intarr(ysize)            ; Y values vector. Undersize
    mean_refcol_l = fltarr(ysize)  ; Signals vector on left
    mean_refcol_r = fltarr(ysize)  ; Signals vector on right
    for y=0, ysize-1 do begin
       row[y] = y

       ; Do not use the first reference pixel on each
       ; output. Rockwell's HAWAII-1RG worked best like this. For
       ; now, assume that HAWAII-2RG will do the same.
       mean_refcol_l[y] = mean(input_data[1:3,y])        ; Average left side reference pixels
       mean_refcol_r[y] = mean(input_data[xsize-4:xsize-2,y]) ; Average right side reference pixels

    endfor

    ; Sigma clip to get rid of outliers
    if (verbose) then $
       print, 'Sigma clipping reference columns to remove outliers'
    sigclip, mean_refcol_l, threshold=3.0, maxiter=5
    sigclip, mean_refcol_r, threshold=3.0, maxiter=5
    if (verbose) then begin
       print, 'Mean reference pixel on left = ', mean(input_data[0:3,*])
       print, 'Mean reference pixel on right = ', mean(input_data[xsize-4:xsize-2,*])
    endif

    ; Apply reference pixel correction
    ; Create arrays to hold columns of values to be subtracted from the science pixels
    refpixcol_left = mean_refcol_l
    refpixcol_right = mean_refcol_r;
    if (not simple) then begin
       ; Apply Spatial Averaging reference pixel correction using
       ; Savitzky-Golay filter. See the SAVGOL help page for
       ; more information.
       ;
       ; Exclude the first and last 4 rows from convolution with the kernel. These are
       ; reference rows and they need to be handled differently for unknown
       ; physical reasons.
       savgol_filter = savgol(fwindow, fwindow, order, degree)
       refpixcol_left[4:ysize-5] = convol(mean_refcol_l[4:ysize-5], savgol_filter, /EDGE_TRUNCATE)
       refpixcol_right[4:ysize-5] = convol(mean_refcol_r[4:ysize-5], savgol_filter, /EDGE_TRUNCATE)
    endif else begin

       ; Apply Maximal Spatial Averaging reference pixel correction. As before,
       ; exclude outtermost reference pixels in columns.
       ;

       ; The approach will be to create a 1D array containing
       ; reference pixel values. Loop over all pixels and load
       ; the reference pixel array

       case scaName of
         'H2RG-003-4.8mu': $ ; For this device, we want to only use the left columns and top and bottom rows in the 1st quadrant
          begin
              ;read reference pixel regions into arrays, reforming into 1-D so they can all be globbed together
              refrows_bottom = reform(input_data[4:(xsize/4)-1, 0:3], n_elements(input_data[4:(xsize/4)-1, 0:3]))
              refrows_top = reform(input_data[4:(xsize/4)-1, ysize-4:ysize-1], n_elements(input_data[4:(xsize/4)-1, ysize-4:ysize-1]))
              refcols_left = reform(input_data[1:3, 4:ysize-5], n_elements(input_data[1:3, 4:ysize-5]))

              ;add all the reference pixels into one array
              allrefpix = [refrows_bottom, refrows_top, refcols_left]
          end
         ELSE: $  ; the default option for all H2RG devices is to just take the average of all reference pixels
          begin
              ;read reference pixel regions into arrays, reforming into 1-D so they can all be globbed together
              refrows_bottom = reform(input_data[4:xsize-5, 0:3], n_elements(input_data[4:xsize-5, 0:3]))
              refrows_top = reform(input_data[4:xsize-5, ysize-4:ysize-1], n_elements(input_data[4:xsize-5, ysize-4:ysize-1]))
              refcols_left = reform(input_data[1:3, 4:ysize-5], n_elements(input_data[1:3, 4:ysize-5]))
              refcols_right = reform(input_data[xsize-4:xsize-2, 4:ysize-5], n_elements(input_data[xsize-4:xsize-2, 4:ysize-5]))

              ;add all the reference pixels into one array
;              allrefpix = [refrows_bottom, refrows_top, refcols_left, refcols_right]
         ; as a test, exclude the bottom row this seems to be bad for the first read
           allrefpix = [refrows_top, refcols_left, refcols_right]
          end
       endcase

       ; Calculate mean of all reference pixels using
       ; djs_iterstat and 3 sigma clipping (the default)
       djs_iterstat, allrefpix, mean=mean_refpix

       ; Both reference columns are the same for Maximal Spatial Averaging
       refpixcol_left[*] = mean_refpix
       refpixcol_right = refpixcol_left

    endelse

    ; Make plots if requested

    if (plot) then begin

       set_device = 'WIN'
 ;      window, retain=2; Necessary on some x displays
       !P.MULTI = [0,1,2]   ; Plot 1 plot over the other

       ; First plot left hand reference pixels
       plot, row, mean_refcol_l, $
        title = 'Mean Reference Column [1:3,*]', $
        xtitle = 'Row Number', $
        ytitle = 'Signal (ADU)', $
        xrange = [0,xsize-1], $
        xstyle = 1

       ; Draw a smooth curve through pixels
       oplot, row, refpixcol_left, thick=3

       ; Now plot right hand reference pixels
       plot, row, mean_refcol_r, $
        title = 'Mean Reference Column [2044:2046,*]', $
        xtitle = 'Row Number', $
        ytitle = 'Signal (ADU)', $
        xrange = [0,xsize-1], $
        xstyle = 1

       ; Draw a smooth curve through points.
       oplot, row, refpixcol_right, thick=3

       !P.MULTI = 0 ; Return to default plotting mode

    endif
;    print, "avg of left ref pixels ",avg(refpixcol_left)
    ;Perform reference pixel correction
;     if (not simple) then begin
        for x=0, (xsize/2)-1 do begin
            output_data[x,*] = input_data[x,*] - refpixcol_left
            output_data[(xsize-1)-x,*] = input_data[(xsize-1)-x,*] - refpixcol_right
        endfor
;     endif else begin
;       output_data = input_data - mean_refpix
;     endelse

end