; NAME:
;			crosstalk_reduce.pro (v1.4)
;
; PURPOSE:
;			Generates Fowler-16 subtracted planes for a list of input images and writes them out to FITS files.  A list
;			of the images produced is generated--this list can be used as the input for crosstalk_analysi.
;
; CALLING SEQUENCE:
;			crosstalk_reduce, 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 'dark.lst'
;			outdir		the name of the directory where output plots are written; defaults to testdir+'Crosstalk'
;			outfile		name of output file list; defaults to 'crosstalkfiles.txt'
;
;
; EXAMPLE:
;			Create output files in a "crosstalk" subdirectory of the current directory:
;			crosstalk_reduce, '.'
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:  Ernie Morse, May 21, 2003
;
;		Ernie Morse, IDTL, May 23, 2003 (v1.1)
;		Procedure now checks each read for cosmic ray hits and writes a cosmic ray mask out to a FITS file
;
;		Ernie Morse, IDTL, June 3, 2003 (v1.2)
;		Added more comments to improve readability
;
;		Ernie Morse, IDTL, June 4, 2003 (v1.3)
;		Changed initial mask so that it contains hit intensity values for cosmic ray hits,
;		instead of just a 1
;
;		Ernie Morse, IDTL, June 5, 2003 (v1.4)
;		Added rejection of hits found in 1st 5 reads, hits that are not brightest hits in neighborhood and hits near image border
;		Added creation of location mask (tells which read # cosmic ray hit struck for each pixel)
;
;		Don Figer, IDTL, October 16, 2004 (v1.4)
;		Replaced pathdelim variable with path_sep(). Added line to make output directory.
;

Pro crosstalk_reduce, testdir, infile=infile, outfile=outfile, outdir=outdir

	; 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 path_sep()) then testdir = testdir + path_sep()

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

	; 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 path_sep()) then outdir = outdir + path_sep()
	file_mkdir,outdir

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

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

	; 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)'

	; loop over the filenames, creating an output Fowler-16 frame for each
	for filecount=0,n_elements(filenames)-1 do begin
		; read in the header for the image in order to get the frame size
		fits_read, testdir+filenames[filecount], 0, imghdr, /header_only
		xsize = long(sxpar(imghdr, 'naxis1'))
		ysize = long(sxpar(imghdr, 'naxis2'))
		numreads = long(sxpar(imghdr, 'naxis3'))
		framesize = xsize * ysize

		; create arrays to hold 3-sigma events, 10-sigma events and the final cosmic ray mask
		; also create a frame to hold the intensity of the hits
		crmask3 = intarr(xsize, ysize)
		crmask10 = intarr(xsize,ysize)
		intensitymask = lonarr(xsize, ysize)
		timemask = intarr(xsize, ysize)

		; create arrays to hold the average of the first and last 16 reads of an input image
		group16_1 = lonarr(framesize)
		group16_2 = lonarr(framesize)


		; now that framesize is know, the frames of the image cube can be read
		fits_read, testdir+filenames[filecount], group16_1, first=0, last=framesize-1
		group16_1 = long(group16_1)
;		fits_read, testdir+filenames[filecount], previousFrame, first=framesize, last=framesize*2-1
;		previousFrame = long(previousFrame)
;		previousFrame = reform(previousFrame, xsize, ysize, /overwrite)
;	  	group16_1 = group16_1 + previousFrame
		previousFrame = reform(group16_1, xsize, ysize)

		for readcount = 1, numreads-1 do begin
			print, readcount
			fits_read, testdir+filenames[filecount], currentFrame, first=framesize*readcount, last=framesize*(readcount+1) - 1
			currentFrame = long(currentFrame)

			; if the current frame is one of the first or last 16, add it to the group array
			if (readcount lt 16) then begin
				group16_1 = group16_1 + currentFrame
			endif

			if (readcount ge (numreads - 16)) then begin
				group16_2 = group16_2 + currentFrame
			endif

			currentFrame = reform(currentFrame, xsize, ysize, /overwrite)

			diffFrame = currentFrame - previousFrame

			sdev = stddev(diffFrame[4:xsize-5, 4:ysize-5])
			med = median(diffFrame[4:xsize-5, 4:ysize-5])

			; keep track of all pixels in current frame - previous frame where the pixel value
			; is > 3 or > 10 sigma of all of the pixels in the difference frame
			outliers = where(diffFrame gt (med + sdev*3))
			crmask3[outliers] = crmask3[outliers] + 1
			outliers = where(diffFrame gt (med + sdev*10))
			crmask10[outliers] = crmask10[outliers]+1

			; write the intensity of the 10-sigma hits into the mask
			intensitymask[outliers] = diffFrame[outliers]
			timemask[outliers] = readcount

			; set value of previous frame for next time through the loop
			previousFrame = currentFrame
		endfor

		; divide by 16 to turn the group arrays into average (they are sum arrays up until this point)
		group16_1 = group16_1 / 16.0
		group16_2 = group16_2 / 16.0

		; make a CDS-subtracted frame for the output
		CDSframe = group16_2 - group16_1
		CDSframe = reform(CDSframe, xsize, ysize, /overwrite)

		; make the initial cosmic ray mask.  Cosmic ray hits are assumed to be those pixels
		; that experience one 10-sigma event with no other events of 3-sigma or greater
		; Also, events in the first 5 reads are rejected because hot pixels that reach saturation
		; within this time frame can be incorrectly identified as hits
		maskpix = where(crmask3 eq 1 and crmask10 eq 1 and timemask gt 4, complement = nonhits)

		; mask already has the intensity of all 10-sigma events, so just set the rest of it to 0
		intensitymask[nonhits] = 0
		timemask[nonhits] = 0

		; do another cut to eliminate cosmic ray hits around the border of the image
		intensitymask[0:4, *] = 0
		intensitymask[xsize-5:*, *] = 0
		intensitymask[5:xsize-6, 0:4] = 0
		intensitymask[5:xsize-6, ysize-5:*] = 0

		timemask[0:4, *] = 0
		timemask[xsize-5:*, *] = 0
		timemask[5:xsize-6, 0:4] = 0
		timemask[5:xsize-6, ysize-5:*] = 0

		; now eliminate any hits that aren't the maximum hit in their neighborhood
		hitpix = where(intensitymask gt 0)

		for pixcount=0, n_elements(hitpix)-1 do begin
			hit = hitpix[pixcount]
			hit_x = hit mod xsize
			hit_y = hit / xsize
			if (intensitymask[hit] ne max(intensitymask[hit_x-1:hit_x+1, hit_y-1:hit_y+1])) then begin
				intensitymask[hit] = 0
				timemask[hit] = 0
			endif
		endfor


		; create output file names
		outFITSname = strmid(filenames[filecount], 0, strpos(filenames[filecount], '.fits')) + '_CDS16.fits'
		outMaskname = strmid(filenames[filecount], 0, strpos(filenames[filecount], '.fits')) + '_CRmask_initial.fits'
		timeMaskname = strmid(filenames[filecount], 0, strpos(filenames[filecount], '.fits')) + '_CRmask_time.fits'

		; write out the output files
		fits_write, outdir+outFITSname, CDSframe, imghdr
		fits_write, outdir+outMaskname, intensitymask, imghdr
		fits_write, outdir+timeMaskname, timemask, imghdr

		; write output FITS file names to output list
		printf, outlun, outFITSname + '   ' + outMaskname
	endfor

	; close the output file list
	free_lun, outlun
End