; NAME:
;			crosstalk_reduce_reag.pro (v1.1)
;
; PURPOSE:
;			Generates output cosmic ray masks and median-subtracted cubes for input REAG crosstalk data
;			Output products are used as input for crosstalk_analysis_reag.pro
;
; CALLING SEQUENCE:
;			crosstalk_reduce_reag, testdir[, infile[, outfile [,outdir]]]
;
; INPUTS:
;			testdir		the name of a directory in which the input data file resides.
;
; KEYWORD PARAMETERS:
;			infile		the name of the input data file.  Defaults to 'reagcubes.txt'
;			outdir		the name of the directory where output plots are written; defaults to testdir+'Crosstalk'
;			outfile		name of output file list; defaults to 'reagmasks.txt'
;
;
; EXAMPLE:
;			Create output files in a "Crosstalk" subdirectory of the current directory:
;			crosstalk_reduce, '.'
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:  Ernie Morse, IDTL, July 23, 2003
;
;		Ernie Morse, IDTL, August 11, 2003 (v1.1)
;		Changed method of detecting cosmic rays
;

Pro crosstalk_reduce_reag, testdir, infile=infile, outfile=outfile, outdir=outdir

	; set directory path delimiter based on operating system in use
	if (!version.os_family eq 'unix') then begin
		pathdelim = '/'
	endif else begin
		pathdelim = '\'
	endelse

	; check to see if input directory name ends in path delimiter--if not, add it to the end
	if (strmid(testdir, strlen(testdir)-1, 1) ne pathdelim) then testdir = testdir + pathdelim

	; set outdir to testdir + 'Crosstalk' if not set by user
	if (not keyword_set(outdir)) then outdir = testdir + 'Crosstalk' + pathdelim

	; check to see if outdirectory name ends in path delimiter--if not, add it to the end
	if (strmid(outdir, strlen(outdir)-1, 1) ne pathdelim) then outdir = outdir + pathdelim

	; set outfile to a default value if it is not set by the calling context
	if (not keyword_set(outfile)) then outfile='reagmasks.txt'

	; set infile to a default value if it is not set by the calling context
	if (not keyword_set(infile)) then infile='reagcubes.txt'

	; open the output file list
	openw, outlun, outdir+outfile, /get_lun

	; read the file names from the input file
	readcol, testdir+infile, filenames, format='(A)'

    for fc=0,n_elements(filenames)-1 do begin
    	; read REAG data cube
        fits_read, testdir+filenames[fc], cube, imghd

        ; get cube dimensions
        cubedim = size(cube)
        xsize = cubedim[1]
        ysize = cubedim[2]
        zsize = cubedim[3]

		; construct a data cube to hold output array
        outcube = fltarr(xsize, ysize, zsize)

        ; begin constructing cosmic ray mask
        thresh = 3.0
		mask = intarr(xsize,ysize,zsize)

		; for each pixel of a detector readout frame, find the median and standard deviation
		; down through the cube of frames
        medimg = median(cube, dimension=3)
        sdimg = fltarr(xsize, ysize)
        for j=0,ysize-1 do begin
            for i=0,xsize-1 do begin
                stack = cube[i,j,*]
                sdimg[i,j] = stddev(stack[where(finite(stack) eq 1)])
            endfor
        endfor

		; find size of one readout frame
        framesize = long(xsize) * ysize

		; indicate position of all pixels with value greater than thresh standard deviations
		; above the median
        for k=0,zsize-1 do begin
        	; construct a median subtracted frame for later output
            outcube[*,*,k] = cube[*,*,k] - median(cube[*,*,k])
            zscore = (cube[*,*,k] - medimg) / sdimg
            hitpix = where(zscore gt thresh)
            if (hitpix[0] ne -1) then begin
                hitpix = hitpix + k*framesize
                mask[hitpix] = 1
            endif
        endfor

		; mask is now a mask of cosmic ray hits in the cube.
		; a 2-column/row border should be masked out

		mask[0:1, *, *] = 0
		mask[xsize-2:xsize-1, *, *] = 0
		mask[*, 0:1, *] = 0
		mask[*, ysize-2:ysize-1, *] = 0

		; construct the file names for the output mask and medaian-subtracted cube and
		; write to output file
		outMaskname = strmid(filenames[fc], 0, strpos(filenames[fc], '.fit')) + $
		    '_CRmask_initial.fits'
		outCubename = strmid(filenames[fc], 0, strpos(filenames[fc], '.fit')) + $
		    '_medsub.fits'
		printf, outlun, outCubename + '  ' + outMaskname



		; write the mask out as a FITS file
		fits_write, outdir + outMaskname, mask

		; remove the bscale keyword from the header and write out the median-subtracted
		; version of the input image cube
		sxdelpar, imghd, 'bscale'
		fits_write, outdir + outCubename, outcube, imghd

    endfor

	; close the output file
	free_lun, outlun

End