; NAME: ; readnoise_reduce (v2.9) ; ; PURPOSE: ; This program reads input image filenames from a file, calculates readnoise ; on pairs of images, and writes the results to a file. ; ; CALLING SEQUENCE: ; IDL> readnoise_reduce, testdir, [outdir, infile, refcor, plottitle, region, ; minexp, novar, varonly, tnint] ; ; INPUTS: ; testdir - string containing name of directory holding list of input ; files and the files themselves ; ; KEYWORD PARAMETERS: ; outdir - string containing name of directory to which output should be ; written (default value is testdir + Readnoise_Results) ; infile - input file name. File must be in directory specified by ; testdir. Default name is 'filelist.txt' ; noref - set to disable reference pixel correction ; ma_refcor - set for maximal averaging reference pixel correction ; plottitle - string containing title for plot. This string is appended to ; 'Readnoise test for ' by the rnplot procedure. Default value ; is the detector name read from the FITS header. ; ; region - a 4-vector of integers specifying the region of a 2D image ; plane over which readnoise is calculated. For an array ; section A[x1:x2, y1:y2], region should be specified as ; [x1, x2, y1, y2]. Default value is assigned based on detector ; type. ; minexp - set if read noise in continuous reads is desired (default is ; total noise in 1000-seconds) ; novar - set if variance vs. time should not be calculated (can be ; used with default or minexp set) ; varonly - set if only variance vs. time should be calculated (value of ; minexp ignored if this keyword is set) ; tnint - exposure time to use for total noise option--default is ; 1000 seconds. ; ; NOTE: The default behavior is to apply spatial averaging reference pixel ; correction. If noref is set, the status of ma_refcor is ignored. ; ; EXAMPLE ; ; IDL> rndir = '\\rabbit\rockwell\SB-304-003\warm1\readnoise.15Nov02' ; IDL> resultsdir = rndir + '\Results' ; IDL> readnoise_reduce, rndir, outdir = resultsdir, /refcor, ; region = [4, 2051, 0, 2047] ; ; REFERENCE: ; ; MODIFICATION HISTORY: ; Written by: ; Mike Telewicz, IDTL ; Modified by: ; Ernie Morse, IDTL, November 15, 2002 ; Added header ; Added call to procedure to create PS and GIF plots ; Updated reference pixel correction to call refsub by B.J. Rauscher ; Eliminated sigma clipping and good pixel mask from input (no longer ; used by readnoise.pro) ; Ernie Morse, IDTL, February 4, 2003 ; Changed to reading input FITS images one plane at a time in order ; to allow for handling of very large files. Commented out call to ; mkdiffimg and instead process CDS and Fowler-sampled images within ; the body of this procedure. ; Ernie Morse, IDTL, March 27, 2003 (v2.0) ; Changed to allow processing on long dark ramps ; Ernie Morse & Don Figer, IDTL, April 10, 2003 (v2.1) ; Changed to use readcol to read all images at once instead of ; reading one-by-one in loop ; Changed name to readnoise_reduce.pro ; Added flush command to force buffered output to file immediately ; Ernie Morse, IDTL, April 14, 2003 (v2.2) ; Changed calls to refsub--must now pass sca name instead of sca type ; Ernie Morse, IDTL, April 16, 2003 (v2.3) ; Fixed problem of railed pixels by marking them as NaN before feeding ; to readnoise.pro ; Bernie Rauscher, IDTL, April 23, 2003 (v2.4) ; Added minexp keyword. ; Ernie Morse, IDTL, April 21 - 27, 2003 (v2.5) ; Added a whole bunch of stuff to support calculating variance vs. ; time ; New keywords novar and varonly added ; Ernie Morse, IDTL, May 1, 2003 (v2.6) ; Corrected some bugs introduced in the change to v2.5 ; Ernie Morse, IDTL, May 5, 2003 (kept version at v2.6) ; Added if-then statement to prevent railedpix from being applied as ; an array index if no railed pixels were present (this came up when ; attempting to do a shorted resistor read noise test.) ; Ernie Morse, IDTL, May 6, 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 + 'ReadNoise_Results' ; Ernie Morse, IDTL, February 5, 2004 (v2.8) ; Edited to ensure compliance with IDTL coding standards ; Sito Balleza, IDTL, March 20, 2004 ; Added SIPIN in detector names. ; Ernie Morse, IDTL, May 18, 2004 (v2.9) ; Added total noise interval keyword (tnint) to define the exposure time ; to use for the total noise test. ; pro readnoise_reduce, testdir, outdir = outdir, infile = infile, $ ma_refcor = ma_refcor, noref = noref, plottitle = plottitle, $ region = region, minexp = minexp, novar = novar, varonly = varonly, $ tnint = tnint ; specify file path delimiter based on operating system 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 begin testdir = testdir + pathdelim endif ; set outdir to testdir if not set by user if (keyword_set(outdir) eq 0) then begin outdir = testdir + 'ReadNoise_Results' + pathdelim endif if (not keyword_set(infile)) then begin infile = 'filelist.txt' 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 ; open files for output if (not keyword_set(varonly)) then begin if (keyword_set(minexp)) then begin datfile = outdir + 'rdnoisevals.txt' plotfile = outdir + 'rdnoiseplot' endif else begin datfile = outdir + 'noisevals.txt' plotfile = outdir + 'noiseplot' endelse endif if (not(keyword_set(novar))) then begin varfile = outdir + 'variancevals.txt' varplot = outdir + 'varianceplot' endif if (not keyword_set(tnint)) then begin tnint = 1000.0 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 ; create non-empty string needed as input to refsub version parameter ; version # returned into ref_ver by refsub 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 ; open input file and read all filenames into an array readcol, testdir + infile, imgfiles, format = '(A)' ; loop through the images in the file list, performing read noise ; on difference images ; paircount keeps track of # of pairs of "like" images processed paircount = 0 ; readsum holds a running sum of the calculated readnoise-- it is divided by ; paircount later to determine average read noise over all the pairs of ; "like" images processed readsum = 0 ; ; for imgcount = 0, n_elements(imgfiles) - 2 do begin filename1 = imgfiles[imgcount] filename2 = imgfiles[imgcount + 1] ; Determine the readmode of both our images (read header only) and get ; some other needed values from the image headers. ; first get the info for the current image fits_read, testdir + filename1, image1, header1, /header_only readmode1 = sxpar(header1, "readmode") sample1 = sxpar(header1, "numsamrd") numutr1 = sxpar(header1, 'numutr') temp1 = sxpar(header1, 'dettemp') naxis1 = sxpar(header1, 'naxis1') naxis2 = sxpar(header1, 'naxis2') naxis3 = sxpar(header1, 'naxis3') exptime1 = sxpar(header1, 'exp_time') ; get info for the next image in the list fits_read, testdir + filename2, image2, header2, /header_only readmode2 = sxpar(header2, "readmode") sample2 = sxpar(header2, "numsamrd") numutr2 = sxpar(header2, 'numutr') temp2 = sxpar(header2, 'dettemp') ; if temp1 equals 0, that means the dettemp header key was not present ; and the experiment was performed without temperature control. ; Assume a warm test (295 K) in this situation if (temp1 eq 0) then begin temp1 = 295 endif ; if this is the first image overall, get some information common to all ; the images from the header if (imgcount eq 0) then begin ;read detector name from header scaName = sxpar(header1, 'detname') ;parse sca type from detector name scaType = (strsplit(scaName, '-', /extract))[0] ; set default region based on scatype if (not keyword_set(region)) then begin case scaType OF 'H1RG' : region = [4, 1019, 4, 1019] 'H2RG' : region = [4, 2043, 4, 2043] 'SIPIN': region = [4, 2043, 4, 2043] 'SB304': region = [4, 2051, 0, 2047] ELSE : begin print, 'Detector type ' + scatype + ' not supported.' $ + ' Exiting procedure.' return end endcase endif ; get parameters for plot routine if (keyword_set(plottitle) eq 0) then begin plottitle = scaName endif ; sxpar doesn't get the date value from IDTL FITS header, so the ; next two lines find the header element that contains the date and ; extract the text after the ' = ' character plotdate = (header1[where(strpos(header1, 'DATE =') ne -1)])[0] plotdate = getwrd(plotdate, 2) ; create string representation of region, to be used later as a ; plotting label region_str = strtrim(string(region), 2) ; write some information to the header of the output file if (not keyword_set(varonly)) then begin openw, lun2, datfile, /get_lun printf, lun2, plottitle printf, lun2, ref_ver printf, lun2, ref_type printf, lun2, plotdate printf, lun2, region_str, $ format = '("[", A, ":", A, ", ", A, ":", A, "]")' endif ; open the file for the variance data, unless the keyword was set ; to disable it ; print header information if (not keyword_set(novar)) then begin openw, varlun, varfile, /get_lun printf, varlun, plottitle printf, varlun, ref_ver printf, varlun, ref_type printf, varlun, plotdate printf, varlun, region_str, $ format = '("[", A, ":", A, ", ", A, ":", A, "]")' endif endif ; now begin processing all the data ; check each image against the next one in the queue--if they have the ; same readmode, # of reads and detector temperature, ; find read noise on their difference image. Otherwise, go on to the ; next image. if (readmode2 eq readmode1 AND sample1 eq sample2 $ AND numutr1 eq numutr2 $ AND round(temp1) eq round(temp2)) then begin ; if we are just starting with images of a certain readmode, number ; of reads and temperature, initialize some values ; readsum is an array that holds a running sum of read noise values ; numreps determines the number of Fowler pairs that should be ; processed from the pair of images ; backstartval is the # of the read representing the 1st read of the ; 2nd group for Fowler sampling ; pairs is an array that holds the Fowler samples that will be ; processed on a pair of images (i.e, [1, 2, 4, 8] to do Fowler-1, ; Fowler-2, Fowler-4 and Fowler-8) if (paircount eq 0) then begin ; for Fowler exposures if (readmode1 ne 'UTR') then begin numreps = 1 readsum = fltarr(1) variancesum = fltarr(2, numutr1) ; find the exposure time for each read, save in first column ; of variance array variancesum[0, *] = (findgen(numutr1) / (numutr1 - 1)) $ * exptime1 backstartval = sample1 pairs = [sample1] utrflag = 0 ; for UTR exposures, simulated Fowler sampling is used endif else begin variancesum = fltarr(2, numutr1) ; find the exposure time for each read, save in first column ; of variance array variancesum[0, *] = (findgen(numutr1) / (numutr1-1)) $ * exptime1 ; store the index number of the read that is closest to ; tnint seconds in the variable targetread ; this will be the read that is the starting point from ; which the 2nd group of Fowlers will be read for ; tnint-second exposure time processing mindiff = min(abs(variancesum[0, *] - tnint), targetread) ; if the minexp keyword is set, the numpairs is set to 16 ; for now. This value is 1/2 the maximum Fowler-N value ; that will be used (so setting it to 16 gives us Fowler-32). if (keyword_set(minexp)) then begin numpairs = 8 endif else begin ; if the keyword is not set, the number of Fowler ; samples that can be done is determined by ; how many reads there are past the target read. ; However, it is currently set to give us no more than ; Fowler-32 because of the excessive processing time ; required to do any more than that. numpairs = ((numutr1 - targetread) / 2) < 16 endelse ; set up the pairs array pairs = intarr(1) + 1 ; pairmult is a multiplicative value to be used shortly pairmult = 1 ; set up the pairs array to hold the number of Fowler pairs ; that we want to process while (pairmult le numpairs) do begin pairmult = pairmult * 2 pairs = [pairs, pairmult] endwhile ;the pairs array should now contain something like ; [1, 2, 4, 8, 16, 32], which is the Fowlers that we want to ; consider numreps = n_elements(pairs) ; we will process all the Fowlers on a pair of images before ; moving on to the next pair.\ ; readsum will hold a running total of the read noise values ; for each Fowler-N under consideration. readsum = fltarr(numreps) endelse endif ; if we're doing variance, set all exposure times to negative to ; indicate that all reads need to be processed if (not keyword_set(novar)) then begin variancesum[0, *] = -1 * ((findgen(numutr1) / (numutr1-1)) $ * exptime1) endif ; define arrays to hold first and last groups of Fowler data for ; like images if (not keyword_set(varonly)) then begin ; 1st Fowler group of 1st image group1_1 = lonarr(naxis1, naxis2) ; 2nd Fowler group of 1st image group1_2 = lonarr(naxis1, naxis2) ; 1st Fowler group of 2nd image group2_1 = lonarr(naxis1, naxis2) ; 2nd Fowler group of 2nd image group2_2 = lonarr(naxis1, naxis2) endif ; stepval determines how many pixels should be read in order to ; read one frame stepval = long(naxis1) * naxis2 - 1 ; this variable will keep track of how many frames have been read ; into each of the group arrays. this is used to prevent ; re-reading the same data multiple times frames_read = 0 ; loop through the various Fowler modes for a pair of images for fowlercount = 0, numreps-1 do begin ; if minexp is set, we need to set a different value for ; backstartval every time through the fowlercount loop ; otherwise, we always use the same value (the read represented ; by targetread) if (keyword_set(minexp)) then begin backstartval = pairs[fowlercount] endif else begin backstartval = targetread endelse ; read planes for 1st group and store their running total in ; arrays ; set beginning point to beginning of frame after the last read ; frame startwd = (stepval + 1) * frames_read ; if minexp is set, the first Fowler group should be read only ; for the Fowler-1 case. ; Otherwise, we can just move add group 2 to group 1 to get the ; new group 1. if (keyword_set(minexp)) then begin if (fowlercount eq 0) then begin group1end = 1 endif else begin group1end = 0 endelse endif else begin group1end = pairs[fowlercount] - frames_read endelse for i = 1, group1end do begin endwd = startwd + stepval ; read from the current read from both images fits_read, testdir + filename1, image1, header1, $ first = startwd, last = endwd fits_read, testdir + filename2, image2, header2, $ first = startwd, last = endwd ; determine the absolute read # readnum = (i + frames_read) - 1 print, 'Processing pair ' + strtrim(string(paircount), 2) $ + ', read # ' + strtrim(string(readnum), 2) $ + ' in 1st Fowler group' ; find the railed pixels in the first read of the 1st image if (frames_read eq 0) then begin railedpix = where(image1 gt 32760 OR image1 lt -32760) firstread1 = image1 firstread2 = image2 endif else begin if (not(keyword_set(novar))) then begin if (readnum ne 0) then begin cds1 = image1 - firstread1 cds2 = image2 - firstread2 cds1 = reform(cds1, naxis1, naxis2, /overwrite) cds2 = reform(cds2, naxis1, naxis2, /overwrite) ; apply reference pixel correction if requested ; by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, cds1, sca = scaName refsub, cds2, sca = scaName endif else begin refsub, cds1, sca = scaName, /simple refsub, cds2, sca = scaName, /simple endelse endif if (railedpix[0] ne -1) then begin cds1[railedpix] = !values.f_nan cds2[railedpix] = !values.f_nan endif varianceval = readnoise(cds1, cds2, region) variancesum[1, readnum] = $ variancesum[1, readnum] + varianceval ; make the time positive to indicate that this ; read has been processed variancesum[0, readnum] = -1 $ * variancesum[0, readnum] endif endif endelse if (not keyword_set(varonly)) then begin ; add the current reads to the group arrays group1_1 = group1_1 + image1 group2_1 = group2_1 + image2 endif startwd = endwd + 1 endfor ;find start read position for 2nd group startwd = backstartval * (stepval + 1) + (stepval + 1) $ * frames_read ; if minexp is set, we should read the whole N reads into group ; 2--if not, we can just read the ones that aren't already there ; this will be accomplished by setting frames_read to 0 if ; minexp is set. this is done later in the code ; read planes for 2nd group and store their running total in ; arrays for j = 1, pairs[fowlercount] - frames_read do begin endwd = startwd + stepval fits_read, testdir + filename1, image1, header1, $ first = startwd, last = endwd fits_read, testdir + filename2, image2, header1, $ first = startwd, last = endwd ; determine the absolute read # readnum = (j + frames_read + backstartval) - 1 print, 'Processing pair ' + strtrim(string(paircount), 2) $ + ', read # ' + strtrim(string(readnum), 2) $ + ' in 2nd Fowler group' if (not(keyword_set(novar))) then begin cds1 = image1 - firstread1 cds2 = image2 - firstread2 cds1 = reform(cds1, naxis1, naxis2, /overwrite) cds2 = reform(cds2, naxis1, naxis2, /overwrite) ; apply reference pixel correction if requested by the ; user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, cds1, sca = scaName refsub, cds2, sca = scaName endif else begin refsub, cds1, sca = scaName, /simple refsub, cds2, sca = scaName, /simple endelse endif if (railedpix[0] ne -1) then begin cds1[railedpix] = !values.f_nan cds2[railedpix] = !values.f_nan endif varianceval = readnoise(cds1, cds2, region) variancesum[1, readnum] = variancesum[1, readnum] $ + varianceval ; make the time positive to indicate that this read has ; been processed variancesum[0, readnum] = -1 * variancesum[0, readnum] endif if (not keyword_set(varonly)) then begin group1_2 = group1_2 + image1 group2_2 = group2_2 + image2 endif startwd = endwd + 1 endfor if (not keyword_set(varonly)) then begin ; construct Fowler-subtracted image for 1st image image1a = (float(group1_2) - float(group1_1)) $ / (pairs[fowlercount]) ; construct Fowler-subtracted image for 2nd image image2a = (float(group2_2) - float(group2_1)) $ / (pairs[fowlercount]) ; apply reference pixel correction if requested by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, image1a, sca = scaName refsub, image2a, sca = scaName endif else begin refsub, image1a, sca = scaName, /simple refsub, image2a, sca = scaName, /simple endelse endif ; calculate the read noise and print the result ; set the railed pixels to NaN before feeding to readnoise if (railedpix[0] ne -1) then begin image1a[railedpix] = !values.f_nan image2a[railedpix] = !values.f_nan endif readnoiseval = readnoise(image1a, image2a, region) readsum[fowlercount] = readnoiseval + readsum[fowlercount] endif ; update the variable frames_read--it should now be equal to the ; Fowler-N value if minexp is not set ; if minexp is set, it should be 0 and the 2nd group should be ; added to the first if (keyword_set(minexp)) then begin frames_read = 0 group1_1 = group1_1 + group1_2 group2_1 = group2_1 + group2_2 group1_2[*] = 0.0 group2_2[*] = 0.0 endif else begin frames_read = pairs[fowlercount] endelse endfor ; do the variance calculation on all reads that didn't get done ; during the Fowler processing if (not keyword_set(novar)) then begin reads_notdone = where(variancesum[0, *] lt 0.0) for varcount = 0, n_elements(reads_notdone) - 1 do begin framenum = reads_notdone[varcount] print, 'Processing pair ' + strtrim(string(paircount), 2) $ + ', read # ' + strtrim(string(framenum), 2) $ + ' post-Fowler' fits_read, testdir + filename1, image1, header1, $ first = framenum * (stepval + 1), $ last = (framenum + 1)*(stepval + 1) - 1 fits_read, testdir + filename2, image2, header1, $ first = framenum * (stepval + 1), $ last = (framenum + 1) * (stepval + 1) - 1 cds1 = image1 - firstread1 cds2 = image2 - firstread2 cds1 = reform(cds1, naxis1, naxis2, /overwrite) cds2 = reform(cds2, naxis1, naxis2, /overwrite) ; apply reference pixel correction if requested by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, cds1, sca = scaName refsub, cds2, sca = scaName endif else begin refsub, cds1, sca = scaName, /simple refsub, cds2, sca = scaName, /simple endelse endif if (railedpix[0] ne -1) then begin cds1[railedpix] = !values.f_nan cds2[railedpix] = !values.f_nan endif varianceval = readnoise(cds1, cds2, region) variancesum[1, framenum] = variancesum[1, framenum] $ + varianceval ; set time back to positive variancesum[0, framenum] = -1 * variancesum[0, framenum] endfor endif paircount = paircount + 1 ; if the two images aren't the same, then continue to the next after ; printing results endif else begin ; print results to file ; paircount is set to 0 after results are printed, so if it is ; now 0, there are no new results for this temperature if (paircount ne 0) then begin if (not keyword_set(varonly)) then begin for rescount = 0, n_elements(pairs)-1 do begin ; print the average read noise values for all pairs ; processed at this detector temperature; ; this is accomplished by dividing the summed readnoise ; by the number of pairs processed. print, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) printf, lun2, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) endfor ;force buffered output to file flush, lun2 endif if (not keyword_set(novar)) then begin ; the value we want to print is the square of the average ; noise from all the pairs processed at this detector ; temperature variancesum[1, *] = (variancesum[1, *] / paircount) ^ 2 printf, varlun, round(temp1), numutr1 printf, varlun, variancesum printf, varlun flush, varlun endif endif paircount = 0 endelse endfor ; if at end of file, print the final result to output file if (paircount ne 0) then begin if (not keyword_set(varonly)) then begin for rescount = 0, n_elements(pairs)-1 do begin ; print the average read noise values for all pairs processed at ; this detector temperature; ; this is accomplished by dividing the summed readnoise by the ; number of pairs processed. printf, lun2, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) print, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) endfor endif if (not keyword_set(novar)) then begin ; the value we want to print is the square of the average noise from ; all the pairs processed at this detector temperature variancesum[1, *] = (variancesum[1, *] / paircount) ^ 2 printf, varlun, round(temp1), numutr1 printf, varlun, variancesum endif endif if (not keyword_set(varonly)) then free_lun, lun2 if (not keyword_set(novar)) then free_lun, varlun if (not keyword_set(varonly)) then begin ; call plot routine to make plots if (keyword_set(minexp)) then begin ; If minexp keyword is set, make read noise plot rnplot, datfile, plotfile, minexp = minexp endif else begin ; Else make total noise plot (the default for RNPLOT.PRO) rnplot, datfile, plotfile, tnint = tnint endelse endif ; make the variance plots if (not keyword_set(novar)) then begin variance_analysis, outdir endif end