; NAME: ; pro darkcurrent_reduce.pro (v2.15) ; ; PURPOSE: ; This program takes a images listed in filename and produces estimates of ; dark current from them. A per-pixel dark current frame is produced, and ; individual dark current slope plots can optionally be created. ; ; CALLING SEQUENCE: ; darkcurrent_reduce, testdir, [outdir [, infile [, valoutfile $ ; [, noref [, ma_refcor [, region [,imgoutfile] [, groupfactor]]]]]] ; ; INPUTS: ; testdir string containing directory name wherein input images and a ; list of the image file names can be found ; ; KEYWORD PARAMETERS: ; outdir directory where output will be written; defaults to ; testdir + '\SlopePlots' ; infile input file name. File must be in directory specified by ; testdir. Default name is 'dark.lst' ; valoutfile name for output file containing read-by-read mode values. ; Will be written to outdir. Default is 'slopevals.txt'. ; imgoutfile name for output file containing names of slope images. ; Will be written to outdir. Default is 'slopeimgs.txt'. ; region rectangular region of SCA to consider for analysis. Should ; be entered as 4-element 1-d integer array, i.e. ; [x1,x2,y1,y2] to consider pixels in rows x1 to x2 and ; columns y1 to y2. Default is specified as all non-reference ; pixels ; ma_refcor set for maximal averaging reference pixel correction ; noref set to disable reference pixel correction ; gdmask set equal to the name of a good pixel mask. 1=good & 0=bad ; groupfactor sets number of reads placed together in a group for ; processing; defaults to 1 (which means mode of each read ; is determined). ; ; ; NOTE: Default behavior is to do spatial average reference pixel correction. ; If both ma_refcor and noref are set, the noref keyword over-rides the ; ma_refcor keyword and no reference pixel correction will be done. ; ; EXAMPLE ; ; Reduce data from specified directory, using default file names, ; generating slope plots and using maximal averaging reference pixel ; correction, using default region. ; ; IDL> dkdir = '\\rabbit\rockwell3\H2RG-006-5.0mu\cold1\dktest.3Apr03' ; IDL> darkcurrent_reduce, dkdir, /ma_refcor ; ; REFERENCE: ; ; ; MODIFICATION HISTORY: ; Written by: ; Don Figer, IDTL, May 3, 2003 ; ; Modified by: ; Don Figer, IDTL, May 5, 2002 ; Added output switches ; ; Bernie Rauscher, IDTL, August 5, 2002 ; Added reference pixel ; correction, glow correction, and spiffy line-fit plots. ; ; Ernie Morse, IDTL, October 21, 2002 (v1.1) ; Modified code so that choosing bjrmods option does ONLY the bjrmods ; section, and not the rest of the procedure ; Added code to create output file for use by darkjumbo2.pro ; Added readtime and month2days functions for use in creating ; darkjumbo output file ; Increased min and max histogram bin limits--they were too small for ; higher-signal images ; ; Ernie Morse, IDTL, October 31, 2002 (v1.2) ; Modified bjrmods section so that PostScript output files are put in ; outdir, not the current directory. ; ; Ernie Morse, IDTL, Nov. 7, 2002 (v1.2.1) ; Now gets SCA type from header, instead of assuming H1RG ; Calls refsub with SCA type instead of refsub_h1rg directly ; ; Ernie Morse, IDTL, January 7, 2003 (v1.2.2) ; Changed BJRMODS histogram setup so that histogram bins are ; centered around median signal in data instead of around 0. ; In the previous method, the actual distribution curve tended to ; move off the edge of the histogram range as signal increased and ; there were a lot of empty bins with negative values. ; ; Ernie Morse, IDTL, February 10-12, 2003 (v1.3) ; Changed reading of FITS files to frame-at-a-time to allow for ; reading arbitrarily large files to prevent limitations imposed by ; the amount of RAM available. ; Added region keyword to allow specification of custom region. ; Removed hard-coded assumptions that array size is 1024 x 1024 ; (array size now taken from header keywords). ; Removed hard-coded assumptions of analysis region. ; Changed plot parameters (increased line & char size & thickness, ; only print det. temp. to 1 mK accuracy). ; ; Ernie Morse, IDTL, March 5, 2003 (v1.3) ; Added keywords to allow specification of type of reference pixel ; correction ; Added information about analysis region and reference pixel ; correction to plots. ; ; Ernie Morse, IDTL, March 13, 2003 (v1.3) ; Added name of file to bottom margin area of slopefit plots in ; bjrmods section. ; ; Ernie Morse, IDTL, April 4, 2003 (v2.0) ; Started all over from scratch. ; ; Ernie Morse, IDTL, April 8, 2003 (v2.1) ; Cast values read in from headers to int and float data types to ; insure that the slope array ; written out at the end will be a float array in order to decrease ; size of the resulting FITS file. ; ; B.J. Rauscher, IDTL, April 10, 2003 (v2.2) ; Added ability to use good pixel masks ; ; Ernie Morse, IDTL, April 11, 2003 (v2.3) ; Changed calls to refsub--must now pass SCA name instead of SCA ; type ; ; Ernie Morse, IDTL, April 14, 2003 (v2.4) ; Removed the bin representing 0 from the gaussfit to protect against ; spurious zeros caused by railed pixels. ; ; Ernie Morse, IDTL, April 16, 2003 (v2.5) ; Corrected bug found in implementation of zero-bin rejection. ; ; Ernie Morse, IDTL, April 30, 2003 (v2.6) ; Added the rd_time keyword and code to allow the processing of the ; Fowler glow images. ; ; Ernie Morse, IDTL, May 9, 2003 (v2.7) ; Eliminated sa_refcor keyword--spatial average reference pixel ; correction is now default ; Added noref keyword--can be set to explicitly disable reference ; pixel correction ; Made changes in code to reflect new default case (spatial averaging). ; Changed default output directory to testdir+'SlopePlots' ; ; Ernie Morse, IDTL, May 15, 2003 (v2.8) ; Added write of reference pixel correction version number and type ; to header of slope images ; ; Ernie Morse, IDTL, October 10, 2003 (v2.9) ; Added cosmic ray detection to creation of per-pixel images ; ; Ernie Morse, IDTL, October 17, 2003 (v2.10) ; Added groupfactor keyword and modifications to allow grouping ; of reads for processing. Split outfile keyword into valoutfile and ; imgoutfile. Made some changes to format of valoutfile (it now ; must include read times as they can not be calculated in the ; analysis script anymore due to the addition of grouping). ; ; Ernie Morse, IDTL, December 2, 2003 (v2.11) ; Corrected bugs associated with operation when groupfactor = 1. ; Added graceful exit for case when all pixels in input image are ; railed. ; ; Ernie Morse, IDTL, July 15, 2004(v2.12) ; Changed from using binsize = 1.0 to starting out with 128 bins with ; size determined by data distribution. ; Added code to catch math errors (from curvefit) ; Corrected problems that arose when first reads of input image are ; at bottom rail. ; Ernie Morse, IDTL, July 22, 2004 (v2.13) ; Added rowinc and colinc keywords. ; Ernie Morse, IDTL, September 29, 2004 (v2.14) ; Modified error handler so that math errors are distinguished from ; non-math errors. ; Ernie Morse, IDTL, October 8, 2004 (v2.15) ; Made changes to prevent blank lines being written to imgoutfile for ; input images that are railed in every read. Moved opening of ; imgoutfile to end of procedure, conditional on there being at least ; one non-railed image and trimmed imginfo array to only include ; entries with non-null file names. ; ;1234567890123456789012345678901234567890123456789012345678901234567890123456789 Pro darkcurrent_reduce, testdir, outdir=outdir, infile=infile, $ valoutfile=valoutfile, imgoutfile=imgoutfile, ma_refcor=ma_refcor, $ noref=noref, region=region, gdmask = gdmask, groupfactor=groupfactor, $ rowinc = rowinc, colinc = colinc ; set directory path delimiter & end-of-line character based on ; operating system in use (non-unix is assumed to be Windows). if (!version.os_family eq 'unix') then begin pathdelim = '/' endofline = string(10b) endif else begin pathdelim = '\' endofline = string(13b) + string(10b) 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 begin testdir = testdir + pathdelim endif ; set outdir to testdir + 'SlopePlots' if not set by user if (keyword_set(outdir) eq 0) then begin outdir = testdir + 'SlopePlots' + pathdelim endif ; 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 begin outdir = outdir + pathdelim endif ; set defaults for keyword values not provided by the calling context if (not keyword_set(infile)) then begin infile = 'dark.lst' endif if (not keyword_set(valoutfile)) then begin valoutfile = 'slopevals.txt' endif if (not keyword_set(imgoutfile)) then begin imgoutfile = 'slopeimgs.txt' endif if (not keyword_set(groupfactor)) then begin groupfactor = 1 endif ; If reference pixel correction is to be done, query the refsub function for ; the version number, otherwise, just set a flag indicating it was not done if (not keyword_set(noref)) then begin ref_ver = ' ' refsub, 0, sca=' ', version = ref_ver if (not keyword_set(ma_refcor)) then begin ref_type = 'Spatial averaging' endif else begin ref_type = 'Maximal averaging' endelse endif else begin ref_ver = 'Not used' ref_type = 'N/A' endelse if (not keyword_set(rowinc)) then begin rowinc = 1 endif if (not keyword_set(colinc)) then begin colinc = 1 endif ; open input file and read all filenames into an array readcol, testdir + infile, imgfiles, format='(A)' ; make an array of structs to hold filename and detector temperature for ; output slope images imginfo = replicate({name:'', dettemp:''}, n_elements(imgfiles)) ; open output file for slope values openw, lun_val, outdir + valoutfile, /get_lun ; start by reading just the header of the first image in the list to get ; some required information fits_read, testdir + imgfiles[0], 0, firsthdr, /header_only xsize = fix(sxpar(firsthdr, 'naxis1')) ysize = fix(sxpar(firsthdr, 'naxis2')) scaName = sxpar(firsthdr, 'detname') ; define the number of pixels to read in for one group of sample reads numpixrd = long(xsize) * ysize ; extract SCA type from SCA name scaType = strmid(scaName, 0, strpos(scaName, '-')) ; if the user didn't set the region, define default regions based on the ; SCA type if (not keyword_set(region)) then begin case scaType OF 'H1RG' : begin region = [4, 1019, 4, 1019] end 'H2RG' : begin region = [4, 2043, 4, 2043] end 'SB304': begin region = [4, 2051, 0, 2047] end ELSE : begin print, 'Detector type ' + scaType + $ ' not supported. Exiting procedure.' return end endcase endif ; construct a string representation of the region array region_str = strtrim(string(region), 2) region_out = '[' + region_str[0] + ':' + region_str[1] + ':' $ + strtrim(string(colinc), 2) + ', ' + region_str[2] + ':' $ + region_str[3] + ':' + strtrim(string(rowinc), 2) + ']' ; write values common to all the images to the header of the output file printf, lun_val, ref_ver printf, lun_val, ref_type printf, lun_val, region_out ; construct arrays for cosmic ray detection mask10sig = intarr(xsize, ysize) ; loop through images listed in input file and process them one frame at ; a time for imgcount = 0, n_elements(imgfiles) - 1 do begin ; To keep it simple, the standard formula from Numerical ; Recipes will be used to calculate slope per pixel in an incremental ; fashion. ; ; slope = N * S[x_i * y_i] - S[x_i] * S[y_i] / ( N * S[x_i^2] - ; S[x_i]^2 ) ; First initialize all sums that we will need. To keep it ; simple, I'll keep some redundant sums (e.g., the ; information already is available to the program) ; declare these running-sum arrays for each file so they'll be zeroed ; at the start each time Sxy = dblarr(xsize, ysize) Sx = dblarr(xsize, ysize) Sy = dblarr(xsize, ysize) Sxx = dblarr(xsize, ysize) ; read the 1st frame fits_read, testdir + imgfiles[imgcount], 0, imghdr, /header_only ; extract information from header exposuretime = double(sxpar(imghdr, 'exp_time')) numreads = fix(sxpar(imghdr, 'numutr')) samples_per_read = fix(sxpar(imghdr, 'numsamrd')) dettemp = float(sxpar(imghdr, 'dettemp')) ; Evidently, the last few reads in the dark current experiment for H2RG-003-SIPIN ; are corrupted. This can be seen by displaying the corresponding images. In order ; to address this problem, I (Don) added the following line of code that throws ; away the last 10 reads. if (scaname eq 'H2RG-003-SIPIN' and numreads gt 12) then numreads=numreads-10 ; find out how many frame groups there are numframes = numreads * samples_per_read numgroups = ceil(float(numframes) / groupfactor) groupends = indgen(numgroups) * groupfactor + (groupfactor - 1) groupends[numgroups - 1] = groupends[numgroups - 1] < (numframes - 1) groupstart = 0 ; create an array containing the time for each group of reads readtimes = dindgen(groupfactor, numgroups) / (numreads - 1) * $ exposuretime ; if groupfactor > 1 and the last group isn't completely filled ; (i.e., the number of reads in the exposure is not evenly ; divisible by the groupfactor), the readtimes array will contain ; exposure times for reads that aren't present because it assumes that ; all groups contain the same number of reads. Mark these invalid ; exposure times as NaN so they can be ignored in processing. overtimes = where(readtimes gt exposuretime) if (overtimes[0] ne -1) then begin readtimes[overtimes] = !values.f_nan endif ; make an array to hold the exposure times for each group. Each ; element is the mean value of the exposure times of the reads ; contained in the group. grouptimes = fltarr(numgroups) for groupnum = 0, numgroups - 1 do begin grouptimes[groupnum] = mean(readtimes[*, groupnum], /nan) endfor ; read the frame groups in until the 1st non-railed group is encountered groupnum = -1 repeat begin groupnum = groupnum + 1 fits_read, testdir + imgfiles[imgcount], firstgroup, $ first = (groupends[groupnum]) * numpixrd, $ last = (groupends[groupnum + 1]) * numpixrd - 1 ; transform the first group of reads into an average framecount = groupends[groupnum + 1] - groupends[groupnum] firstgroup = reform(float(firstgroup), xsize, ysize, framecount) firstgroup = total(firstgroup, 3) / framecount regionpix = firstgroup[region[0]:region[1]:colinc, $ region[2]:region[3]:rowinc] lowrailpix = where(regionpix lt -32760) hirailpix = where(regionpix gt 32760) lowfraction = n_elements(lowrailpix) / float(n_elements(regionpix)) hifraction = n_elements(hirailpix) / float(n_elements(regionpix)) hi_or_low = hifraction ge lowfraction ? 32767.0 : -32768.0 endrep until ((lowfraction le 0.5 and hifraction le 0.5) $ OR (groupnum eq (numgroups - 2))) ; skip to next image if all reads are railed low if (groupnum eq numgroups - 2) then begin ; make an array to hold the output combination of group time and ; corresponding mode value and raw signal outcombo = fltarr(3, numgroups) outcombo[0, *] = grouptimes outcombo[2, *] = hi_or_low ; write out time and signal values for each read to the file, ; darkcurrent_analyze can be called later to make slope plots. printf, lun_val, testdir+imgfiles[imgcount] + ' ' + $ strtrim(string(numgroups), 2) + ' ' + $ strtrim(string(groupfactor), 2) printf, lun_val, outcombo, format='(F14.3, F14.3, F14.3)' printf, lun_val ; print blank space to separate data for ; different images continue endif ; the first group is the zero-level for the dark current, so adjust ; the grouptimes to be relative to the time of the first group firstgroupnum = groupnum starttime = grouptimes[firstgroupnum] grouptimes = grouptimes - starttime ; find out where the railed pixels are railedpix = where(firstgroup gt 32760 OR firstgroup lt -32760) ; exit the procedure if every pixel in the read is railed if (n_elements(railedpix) eq n_elements(firstgroup)) then begin stop, 'Processing cannot continue because every pixel in ' + $ 'input image is railed at the top or bottom of the A-D range.' endif ; declare array to hold median and mode values theMode = (theMedian = fltarr(numgroups)) ; declare array to hold raw median signal for each read group rawsig = fltarr(numgroups) rawsig[firstgroupnum] = median(firstgroup[region[0]:region[1]:colinc, $ region[2]:region[3]:rowinc]) if (firstgroupnum ne 0) then begin rawsig[0:firstgroupnum - 1] = -32768 endif ; set up array used to hold previous read for cosmic ray detection previousgroup = fltarr(xsize,ysize) mask10sig[*] = 0 ; loop through the reads of an individual input file for groupcount = firstgroupnum + 1, numgroups - 1 do begin ; read all the frames composing one sample group fits_read, testdir + imgfiles[imgcount], scagroup, $ first = (groupends[groupcount - 1] + 1) * numpixrd, $ last = (groupends[groupcount] + 1) * numpixrd - 1 ; transform the group of reads into an average framecount = groupends[groupcount] - groupends[groupcount - 1] scagroup = reform(float(scagroup), xsize, ysize, framecount) scagroup = total(scagroup, 3) / framecount ; record the median signal rawsig[groupcount] = median(scagroup[region[0]:region[1]:colinc, region[2]:region[3]:rowinc]) ; subtract off the first group as bias correction scagroup = scagroup - firstgroup ; if it was requested, perform reference pixel correction on all ; reads if (not keyword_set(noref)) then begin print, 'Applying reference pixel correction to ' + $ imgfiles[imgcount] + ', read # ' + $ strtrim(string(groupcount),2) if (not keyword_set(ma_refcor)) then begin refsub, scagroup, scaName=scaName endif else begin refsub, scagroup, scaName=scaName, /simple endelse endif else begin print, 'No ref. pixel correction for ' + imgfiles[imgcount] + $ ', read #' + strtrim(string(groupcount),2) endelse ; do cosmic ray detection diffFrame = scagroup - previousgroup crsdev = stddev(diffFrame[4:xsize-5, 4:ysize-5]) crmed = median(diffFrame[4:xsize-5, 4:ysize-5]) outliers = where(diffFrame gt (crmed + crsdev*10)) if (outliers[0] ne -1) then begin mask10sig[outliers] = mask10sig[outliers] + 1 endif previousgroup = scagroup ; do the incremental slope calculation operations Sxy = Sxy + grouptimes[groupcount] * scagroup Sx = Sx + grouptimes[groupcount] Sy = Sy + scagroup Sxx = Sxx + grouptimes[groupcount] ^ 2 ; begin the process of finding signal values for individual reads theMedian[groupcount] = $ median(scagroup[region[0]:region[1]:colinc, $ region[2]:region[3]:rowinc]) sdev = stddev(scagroup[region[0]:region[1]:colinc, $ region[2]:region[3]:rowinc]) ; We want to find the mode. The easiest way to do this is to fit a ; gaussian to a histogram. ; set min and max histogram bins, centered around the median of the ; data ;_binsize = 1.0 ;_min = long((theMedian[groupcount] - (sdev * 4))) + _binsize / 2.0 ;_max = long((theMedian[groupcount] + (sdev * 4))) - _binsize / 2.0 _min = long((theMedian[groupcount] - (sdev * 4))) _max = long((theMedian[groupcount] + (sdev * 4))) _nbins = 128 ; Only the non-railed pixels should be considered, so set the ; railed ones to NaN if (railedpix[0] ne -1) then begin scagroup[railedpix] = !values.f_nan endif ; Make a histogram using either all of the data, or only data where ; gdmask says to. math_err_stat = check_math(mask = 208) math_errcount = 0 catch, errval if (errval ne 0) then begin ; if there is a math error which has been thrown manually ; by a call to the message procedure, reduce the number of bins ; and retry if (strpos(!error_state.msg, 'math error') ne -1) then begin _nbins = _nbins / 2 message, /reset math_errcount = math_errcount + 1 endif else begin ; if this is an unexpected non-math error, cancel the ; error handler and retry, allowing the error to be ; handled by the default IDL handler (which should exit the ; procedure). catch, /cancel endelse endif if (not keyword_set(gdmask)) then begin ; If not using the mask, look at all pixels ;theHist = $ ;histogram(scagroup[region[0]:region[1], $ ;region[2]:region[3]], binsize=_binsize, min=_min, $ ;max=_max, /nan) theHist = histogram(scagroup[region[0]:region[1]:colinc, $ region[2]:region[3]:rowinc], nbins = _nbins, $ min = _min, max = _max, /nan, locations = theVals) endif else begin ; If using the mask, look only at masked pixels in ; photo-sensitive region fits_read, gdmask, gdmask_data, gdmask_header, /data_only ;theHist = histogram(scagroup[where(gdmask_data eq 1)], $ ; binsize=_binsize, min=_min, max=_max, /nan) theHist = histogram(scagroup[where(gdmask_data eq 1)], $ nbins = _nbins, min = _min, max = _max, /nan, $ locations = theVals) endelse ; determine the corresponding values for the histogram bins ;theVals = findgen(n_elements(theHist)) * _binsize + _min + $ ; _binsize / 2.0 ; adjust values to center of bins theVals = theVals + ((_max - _min) / float(_nbins - 1)) / 2.0 ; fit the histogram counts and values to a gaussian theResult = gaussfit(theVals, theHist, A) math_err_stat = check_math(mask = 208) if (math_err_stat ne 0 and math_errcount lt 5) then begin message, 'math error' endif ; we want the mode, which is element 1 of the A array returned from ; gaussfit print, math_errcount, a[1] theMode[groupcount] = A[1] endfor ; make an array to hold the output combination of group time and ; corresponding mode value and raw signal outcombo = fltarr(3, n_elements(theMode)) outcombo[0, *] = grouptimes + starttime outcombo[1, *] = theMode outcombo[2, *] = rawsig ; write out time and signal values for each read to the file, ; darkcurrent_analyze can be called later to make slope plots. printf, lun_val, testdir+imgfiles[imgcount] + ' ' + $ strtrim(string(numgroups), 2) + ' ' + $ strtrim(string(groupfactor), 2) printf, lun_val, outcombo, format='(F14.3, F14.3, F14.3)' printf, lun_val ; print blank space to separate data for ; different images ; calculate the final per-pixel slope values outarray = (numgroups * Sxy - Sx * Sy) / $ (numgroups * Sxx - Sx * Sx) ; change output array to float to make output FITS file smaller outarray = float(outarray) ; set the detected cosmic rays/hot pixels to NaN if (mask10sig[0] ne -1) then begin crpix = where(mask10sig ne 0) if (crpix[0] ne -1) then begin outarray[crpix] = !values.f_nan endif endif ; correct values in the original header to indicate 1 frame of float ; values sxaddpar, imghdr, 'naxis3', 1 sxaddpar, imghdr, 'bitpix', -32 sxaddpar, imghdr, 'refver', ref_ver sxaddpar, imghdr, 'reftype', ref_type ; write out slope array with filename equal to original filename ; with '_slope' inserted prior to '.fits' nameparts = strsplit(imgfiles[imgcount], '.', /extract) outname = strjoin(nameparts[0:n_elements(nameparts)-2], '.') + $ '_slope.fits' fits_write, outdir + outname, outarray, imghdr ; enter the output name and detector temperature (rounded to nearest ; 0.1 K) into imginfo array imginfo[imgcount].name = outname imginfo[imgcount].dettemp = string(dettemp, format='(F5.1)') endfor ; get rid of null file names (representing images that are railed in every read) validnames = where(imginfo.name ne '', numvalid) if (numvalid gt 0) then begin imginfo = imginfo[validnames] ; open output file for names of slope images openw, lun_img, outdir + imgoutfile, /get_lun ; sort the imginfo struct array by the detector temperature field imginfo = imginfo[sort(imginfo.dettemp)] ; find the lengths of all the filenames and set the length of the field for ; all names to the maximum value + 5 characters (for output formatting) namelen = strlen(imginfo.name) fieldlen = max(namelen) + 5 for i=0, n_elements(imginfo) - 1 do begin imginfo[i].name = imginfo[i].name + string(make_array(fieldlen - $ strlen(imginfo[i].name), /byte, val=32b)) endfor ; add EOL character(s) to output dettemp so that output will not all be on ; the same line imginfo.dettemp = imginfo.dettemp + endofline ; write the imginfo struct array out to the image list file writeu, lun_img, imginfo.name + endofline free_lun, lun_img endif ; close the output files free_lun, lun_val End