; 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.
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 H4RG 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]'

    ; 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=indgen(ysize)
    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
       mean_refcol_l[y] = float(mean(input_data[0:3,y]))        ; Average left side reference pixels
       mean_refcol_r[y] = float(mean(input_data[xsize-4:xsize-1,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-1,*])
    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 H4RG 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 [4091:4095,*]', $
        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,*] = float(input_data[x,*]) - float(refpixcol_left)
;            output_data[(xsize-1)-x,*] = float(input_data[(xsize-1)-x,*]) - float(refpixcol_right)
;        endfor
;     endif else begin
;       output_data = input_data - mean_refpix
;     endelse

	; the reference data are taken to be the average of the right and left reference averages
	refpixcol=(refpixcol_right+refpixcol_left)/2.
	for x=0,xsize-1 do begin
		output_data(x,*)=float(input_data[x,*]) - float(refpixcol)
	endfor

end