;+ ; NAME: mktwiflat ; PURPOSE: ; makes a flatfield from a series of images of the twilight sky ; using a iterative fitting technique for each pixels to measure ; the relative flatfield and to account for any constant (i.e. thermal) ; component. to be used for making infrared twiflats. ; ; ; INPUTS: ; KEYWORDS: ; infile name of a text file containing a list of files. ; if present, this overrides istart and iend ; datadir file directory ; sort_em sort 'em! Necessary if files are not in order. ; noflag don't try to flag negative pixels with BADVAL. ; ; OUTPUTS: ; ; HISTORY: ; - ; flags pixels with <0.001 as BADVAL too ; 11/12/98 MCL ; ; small big fix: did not properly flag pixels with ; too few fittable pixels as BADVAL ; now prints (c0,c1) on plots ; 11/23/98 MCL ; ; displays stack of fitting masks (can see if there were stars) ; 01/27/99 MCL ; ; add flag for IRCAL ; ** placed under RCS, version 1.0 ** ; 08/04/99 MCL ; ; plots the residuals to the fit in a 2nd window ; only overplots the final fitted line (instead of each iteration, ; which was basically indistinguishable) ; only prints info for pixels which go through >= 2 iterations ; (used to be >= 1 iteration) ; 08/04/99 MCL ; ; 2002-03-30 23:24:25 Converted to .pro by Marshall Perrin ; only works with IRCAL files for now. ; 2002-04-11 15:33:59 added images keyword to allow passing ; in images directly. Also head0 for header. ; 2002-07-04 15:46:49. Converted to use mloadims ; added "auto" keyword to allow fully automatic execution. ; 2002-12-03 12:21:13. Updated to use fitsloader. ; 2004-04-16 MDP. Changed bad mask code to use 4 sigma for hot pixels ; instead of 3 sigma, which tended to mark too many pixels as bad. ; 2004-05-24 MDP. Added /MPFIT and GAIN= keywords. ; ;---------------------------------------- PRO mktwiflat,istart,iend,infile=infile,biasfile=biasfile,$ fdir=fdir,datadir=datadir,outdir=outdir,sort_em=sort_em,noflag=noflag,$ polmask=polmask,outname=outnam,images=images,header=header,auto=auto,$ caldir=caldir,instrument=instrument,hdrs=hdrs,$ mpfit=mpfit,gain=gain if ~(keyword_set(outdir)) then outdir = './' if keyword_set(fdir) then datadir = fdir ; back-compatibility if ~(keyword_set(datadir)) then datadir = './data' if ~(keyword_set(outnam)) then outnam = 'twiflat' if ~(keyword_set(gain)) then gain = 10.0 ; IRCAL gain e-/DN checkdir,outdir checkdir,datadir if ~(keyword_set(instrument)) then begin ctio = 0 ; if CTIO files gemini = 0 ; if Gemini files lick = 0 ; if Lick files ircal = 1 ; if IRCAL files aobir=0 ; if CFHT AOBIR files nirc2=0 ; if Keck NIRC2 files if (ircal gt 0) then instr="ircal" else instr="kermit" endif disp = 1 ; display the fit (& outliers)? DSTEP = 400 ; plot the fit for every DSTEP'th pixel NITER = 3 ; makes little difference NSIGREJ = 3.5 ; don't mess with this! MINTWI = 3 ; min number of frames needed for fit ;---------------------------------------- BADVAL = -1e6 ;------------------------------------------------------------ i0 = ISTART i1 = IEND ; get images print, '=== loading images ===' if not keyword_set(images) then begin fitsloader,transpose([[i0],[i1]]),instrument=instrument,datadir=datadir,caldir=caldir,$ twi,head0,/firstheader,/dobias,imlist=imlist endif else begin twi = images head0= hdrs[*,0] sxaddhist,"MKTWIFLAT: data passed via images/hdrs arguments.",head0 sz = size(imgs0) endelse ;------------------------------------------------------------ sz=size(twi) imgsize=sz[1] lincoeff = fltarr(imgsize, imgsize) ; history string sxaddhist, 'MKTWIFLAT.IDL: '+systime()+ userid(), head0 if ~(keyword_set(auto)) then stop ; subtract bias ; This is now done by fitsloader, above twi_nb = temporary(twi) ;if keyword_set(biasfile) then begin ; print, '=== subtracting bias ===' ; skysub, twi, bias, twi_nb, /silent ; twi = 0 ; sxaddhist, ' SKYSUB: subtracting bias frame "'+biasfile+'"', head0 ;endif else begin twi_nb=temporary(twi) ;endelse ; selected subset of images to use sz = size(twi_nb) ncols = sz(1) nrows = sz(2) ntwi = sz(3) print, '=== image statistics (bias-subtracted) ===' stats, twi_nb print, 'there are ', strc(ntwi), ' twi frames' while getyn('use only a subset of the images?', 0,noquery=auto) do begin sz = size(twi_nb) nn = sz(3) a0 = -1 & a1 = -1 repeat begin read, 'enter starting image: (0-',+strc(nn-1)+') ', a0 read, 'enter ending image: (0-',+strc(nn-1)+') ', a1 endrep until ((a0 ge 0) and (a0 lt nn) and $ (a1 ge 0) and (a1 lt nn) and (a1 gt a0)) twi_nb = twi_nb(*, *, a0:a1) print, '=== examine image stats (w/o bias) ===' stats, twi_nb i0 = i0 + a0 i1 = i0 + (a1-a0) sxaddhist, ' selected images '+strc(i0)+' to '+strc(i1), head0 ntwi = a1-a0+1 print, 'there are ', strc(ntwi), ' twi frames' endwhile ; apply linearity correction ; (use /nocorrect flag if no lin coeffs - will still check for saturation) if getyn('apply linearity correction?', (avg(lincoeff) ne 0),noquery=auto) then begin lincorr, twi_nb, twi_lin, lincoeff, satmap if total(lincoeff) ne 0 then $ sxaddhist, ' LINCORR: applied linearity correction', head0 $ else $ sxaddhist, ' LINCORR: no linearity correction', head0 endif else begin if not(keyword_set(noflag)) then $ lincorr, twi_nb, twi_lin, fltarr(imgsize, imgsize), satmap, /nocorrect $ else twi_lin = temporary(twi_nb) sxaddhist, ' LINCORR: no linearity correction', head0 endelse twi_nb = 0 SPT: ; extract image medians ; make sure images are sorted in order of increasing median. ; this is necessary because if we combine dawn and dusk flats then ; the linear fit will go bonkers unless they are reordered. twi_meds = meds(twi_lin) if keyword_set(sort_em) then begin s = reverse(sort(twi_meds)) twi_meds=twi_meds[s] twi_lin=twi_lin[*,*,s] endif ;deal with multiple polarization regions. We have to do this -after- ; the step above so we can sort the images first. if (keyword_set(polmask)) then $ twi_medarr = pol_meds(twi_lin,polmask,/join) ; initialize output arrays twi_flat = fltarr(sz(1), sz(2)) twi_const = fltarr(sz(1), sz(2)) twi_masks = fltarr(sz(1), sz(2), sz(3))+1B ; do pixel by pixel fit if (disp) then begin lincolr, /silent win, 0, xs = 400, ys = 600 endif for i = 0, ncols-1 do begin ;print, form = '($,A," ")', strc(i) for j = 0, nrows-1 do begin ; initialize cc = [BADVAL, BADVAL] !p.multi = [0, 1, 2] ; extract counts cts = reform(twi_lin(i, j, *)) w = where(cts ne BADVAL and finite(cts), nw) if (nw eq 0) or (nw lt MINTWI) then $ ni = NITER $ else begin cts = cts(w) if keyword_set(polmask) then tmeds = twi_medarr(i,j,w) $ else tmeds = twi_meds(w) cts0 = cts ni = 0 ; plot counts versus medians if (disp) and (((i*sz(1)+j) mod DSTEP) eq 0) then begin if keyword_set(polmask) then tplot = twi_medarr(i,j,*) $ else tplot = twi_meds ; don't try to plot if there's no valid data here... if total(finite(cts0)) le 2 then ctsplot=tplot*0 else ctsplot=cts plot, tplot, ctsplot, ps = 4, sym = 2, $ xtit = 'median counts', ytit = 'pixel counts', $ tit = 'pixel '+printcoo(i, j) endif endelse ; now do iterative fitting ncc = n_elements(cts) while (ni lt NITER) do begin ; do robust fit to extract gain if keyword_set(mpfit) then begin skynoise = sqrt(cts/gain) cc = mpfitexpr('P[0]+P[1]*X',tmeds,cts,skynoise,[0.,1.0]) endif else begin ; FIXME the following was changed to linfit ; for NIRC2 because ladfit was locking up. cc = ladfit(tmeds, cts) ; trying linfit just fit kicks 2004-05-11 ;cc = linfit(tmeds, cts) endelse ;stop ; calculate deviations from fit and robust sigma ctsfit = cc(0)+cc(1)*tmeds diff = cts-ctsfit iterstat, diff, istat, nsig = NSIGREJ, /silent sigma = istat(3) ; flag outliers ; after each iteration, tolerance of outliers increases w = where(diff le (NSIGREJ+ni)*sigma, nw) if (nw lt ncc) then begin wbad = where(diff gt (NSIGREJ+ni)*sigma) twi_masks(i, j, wbad) = 0 if (disp) and (((i*sz(1)+j) mod DSTEP) eq 0) then begin if total(finite(cts(wbad))) ge 2 then $ oplot, tplot(wbad), cts(wbad), ps = 2, sym = 2, col = ni+2 endif ni = ni+1 ; print info for any pixel which goes through >=2 iterations if (ni ge 2) then begin print, printcoo(i, j), ': ni=', strc(ni), $ ', #out = ', strc(ncc-nw), $ ', wbad = '+commalist(wbad), $ ', Nsig = ', strc(diff(wbad)/sigma) endif ;if (ni gt 1) or (max(abs(diff(wbad))/sigma) gt 10) then $ ; if not(getyn('continue?')) then stop ncc = nw cts = cts(w) tmeds = tmeds(w) endif else $ ni = NITER endwhile ; if not a bad pixel if (nw eq 0) or (nw lt MINTWI) then begin ; do nothing endif else begin ; if desired, overplot the final fit ; and the residuals from the fit if (disp) and (((i*sz(1)+j) mod DSTEP) eq 0) then begin oplot, tmeds, ctsfit, col = 3 ;oplot,tmeds,ctsfit = cc(0)+cc(1)*tmeds xyouts, /data, avg(!x.crange), avg(!y.crange)*0.2, $ align = 0.5, chars = 1.5, $ 'c1 = '+strc(cc(1))+'!Cc0 = '+strc(cc(0)) plot, tmeds, cts-ctsfit, psym = 4, sym = 2, $ xrange = !x.crange, xtit = 'median counts', $ ytit = 'residual = (data) - (fit)' if n_elements(wbad) gt 0 then begin ;oplot, twi_meds(wbad), intarr(n_elements(wbad)), $ oplot, tplot(wbad), $ cts0(wbad)-cc(0)+cc(1)*tplot(wbad), $ ps = 2, sym = 2, col = 2 endif plots, !x.crange, [0, 0], line = 5 ; oplot the expected photon noise if keyword_set(ircal) then begin plots,!x.crange,sqrt(!x.crange/10),line=1,/clip ; ircal gain is 10 e-/DN plots,!x.crange,(-1.)*sqrt(!x.crange/10),line=1,/clip ; ircal gain is 10 e-/DN endif legend, /bottom, /right, ['good', 'bad'], $ psym = [4, 2], col = [!p.color, 2], box = 0 endif endelse twi_flat(i, j) = cc(1) twi_const(i, j) = cc(0) ;print,printcoo(i,j) if (disp) and (((i*sz(1)+j) mod DSTEP) eq 0) then $ print, '** plotting fit for ', printcoo(i, j), ' **' endfor endfor sxaddhist, ' doing pixel-by-pixel robust linear fit', head0 sxaddhist, ' max # of iterations = '+strc(NITER), head0 sxaddhist, ' # of sigma for clipping, NSIGREJ = '+strc(NSIGREJ), head0 sxaddhist, ' min # of frames for fitting = '+strc(MINTWI), head0 !p.multi = [0, 1, 1] ; print number of reject pixels per image print, '--------------------------------------------------' for i = 0, ntwi-1 do begin w = where(twi_masks(*, *, i) eq 0, nw) print, 'image '+strc(i)+': # rejected pixels = ', strc(nw) endfor ; normalize and flag bad pixels (<0.001) twi_flatraw = twi_flat twi_flat = normal(twi_flat) w = where(twi_flat le 0.001, nw) if (nw gt 0) then begin twi_flat(w) = !values.f_nan sxaddhist, ' setting <=0.001 pixels to NaN, # = '+strc(nw), head0 print, ' setting <=0.001 pixels to NaN, # = '+strc(nw) endif ;---------------------------------------- ; guess a bad pixel mask ; IRCAL traditionally has about 1.5% of the pixels low due to the two blocked ; corners and other misc defects ; ; Let's throw away everything more than 3 standard deviations down from the mean ; ; grab cold pixels based on slopes ; print, "Now marking bad pixels " ;binsize = 0.001 ;h = histogram(twi_flat > 0,bin=binsize,omin=omin,omax=omax) ;hi = total(h,/cumulative) ;cutoff = binsize*(value_locate(hi,0.02*n_elements(twi_flat))+1) wg = where(twi_flat ne BADVAL and not finite(twi_flat,/NaN)) cutoff = mean(twi_flat[wg]) - 3*stddev(twi_flat[wg]) sz = size(twi_flat) ; convention is 1 = good, 0 = bad badpixmask = bytarr(sz[1],sz[2])+1B cold = where(twi_flat le cutoff,wcold) if wcold gt 0 then badpixmask[cold] = 0B print, "Cold pixel cutoff: slope lt "+strc(cutoff)+", "+strc(wcold)+" pixels rejected." ;grab hot pixels based on constant offset ; there aren't too many of these.. ;binsize = .5 ;h = histogram(twi_const > 0,bin=binsize,omin=omin,omax=omax) ;hi = total(h,/cumulative) ;cutoff = binsize*(value_locate(hi,0.997*n_elements(twi_flat))) ; actually this new code grabs both hot & cold 3 sigma constant offsets wg = where(twi_const ne BADVAL and not finite(twi_const,/NaN)) mconst = mean(twi_const[wg]) badpixmask[where(abs(twi_const-mconst) ge 4*stddev(twi_const[wg]),whot)] = 0B print, "Hot pixel cutoff: const gt "+strc(4*stddev(twi_const[wg]))+", "+strc(whot)+" pixels rejected." ; now just grab badvals wb = where(twi_flat eq BADVAL,wbv) if wbv gt 0 then badpixmask[wb] = 0B print, "bad pixel: slope is BADVAL, "+strc(wbv)+" pixels rejected." wb = where(twi_const eq BADVAL,wbv) if wbv gt 0 then badpixmask[wb] = 0B print, "bad pixel: const is BADVAL, "+strc(wbv)+" pixels rejected." wbad = where(badpixmask eq 0B,badcount) print, "TOTAL BAD PIXELS: "+strc(badcount) ; display results loadct, 15 win, 0 display2, twi_flat, tit = 'final twiflat' win, 1 display2, twi_const, tit = 'constant emission' iterstat, twi_const, istat, /silent print print, '----------------------------------------' print, ' output statistics' print, '----------------------------------------' print, 'constant frame median = ', strc(median(twi_const)) print, 'constant frame (iter)sigma = ', strc(istat(3)) print, 'constant frame sigma = ', $ strc(stdev(twi_const(where(twi_const ne BADVAL)))) print stack, twi_masks, masktot win, 2 display2, /tv, masktot eq 1, tit = 'total of fitting masks' win,3 display2,badpixmask,title="Bad Pixel Mask (dark = bad)" ; save images if desired if getyn('save the flat?', 1,noquery=auto) then begin if keyword_set(aobir) then begin ; this is needed to change the FITS header since the ; input stuff from CFHT has BZERO=32768 but our output here doesn't. print,"Updating BZERO keyword in output file" sxaddpar,head0,"BZERO",0,"offset" endif if not(keyword_set(auto)) then $ if getyn('outdir="'+outdir+'", outnam="'+outnam+'": '+ $ 'do you want to change these?', 0) then begin read, 'enter outdir: ', outdir read, 'enter outnam: ', outnam endif ; just use first image to name flat. if n_elements(i0) gt 1 then i0 = i0[0]; outfile = outdir+outnam+"."+strc(i0)+'.fits' objname = sxpar(head0, 'OBJECT') if filecheck(outfile) then $ writefits, outfile, twi_flat, head0 outfile = outdir+outnam+"."+strc(i0)+'_masks.fits' sxaddpar, head0, 'OBJECT', objname+' - MASKS' ; the masks file is very large and is usually not consulted ; don't save it during automatic mode if not(keyword_set(auto)) then begin if getyn("Save pixel masks for individual frames? (large file) ",0) then $ if filecheck(outfile) then $ writefits, outfile, twi_masks, head0 endif outfile = outdir+outnam+"."+strc(i0)+'_const.fits' sxaddpar, head0, 'OBJECT', objname+' - CONSTANT' if filecheck(outfile) then $ writefits, outfile, twi_const, head0 outfile = outdir+outnam+"."+strc(i0)+"_badmask.fits" sxaddpar, head0, 'OBJECT', objname+' - BAD PIXELS' if filecheck(outfile) then $ writefits, outfile, badpixmask, head0 endif print print, '** MKTWIFLAT.PRO done **' print end