;************** subpixel version, no clicking on stuff necessary ********* pro redircalsub,objn,istart,iend,pclip,$ fdir=fdir,datadir=datadir,flatfile=flatf,badpixf=badpixf,outdir=outdir,subpixel=subpixel,$ subshift=subshift,submethod=submethod,imgs0=imgs0,biasfile=biasfile,$ auto=auto,medsky=medsky,median=median,nan=nan,saturation=saturation,$ scaletime=scaletime,fixpix=fixpix,shiftfile=shiftfile,caldir=caldir,$ satmask=satmask,instrument=instrument,skyfile=skyfile,submedians=submedians,$ verbosename=verbosename,listfile=listfile,$ dewarp=dewarp ;+ ; ; FUNCTION: redircalsub ; PURPOSE: ; Reduction of IRCAL data with subpixel registration and mosaicing. ; NOTES: ; After the standard bias subtraction and flat fielding steps, ; this procedure uses cross-correlations to register any arbitrary dither ; together and combines the images into a final mosaic. ; ; INPUTS: ; objn Object Name ; istart file name index start ; iend file name index end ; pclip percentage of best images to keep after FWHM cut ; ; KEYWORDS: ; datadir= data directory ; flatf= flat field to use ; badpixf= bad pixel mask file ; outdir= output data directory (must already exist) ; /subpixel do subpixel registration ; subshift= subpixel shift method (see mosf.pro) ; submethod= subpixel registration method (see subreg.pro) ; auto Perform all tasks automatically without prompting the user. ; /medsky create median frame for sky and subtract it from all frames. ; /median use median combining, not averaging. ; nan fix nans before registering, the quick way. ; scaletime for fitsloader. scale images of different exposure times ; to this exposure time. ; saturation also for fitsloader. what value (for an unscaled image) is saturated? mask as NaNs ; /satmask use ircal_satmask function to mask saturated pixels ; /verbosename Put the expose time and istart into the output file name ; ; OUTPUTS: ; ; HISTORY: ; ; quick redux of Lick AO data ; 06/11/01 M. Liu ; 06/13/01 mods by JPL to add flatfield ; 2001-07-17 turned into a procedure by mperrin ; 2001-08 subpixel stuff added ; 2002-04-09 dark frame subtraction added by mperrin ; 2002-07-04 nstack parameter removed - not used by the subpixel code ; added auto keyword. ; 2002-08-06 added medsky keyword ; 2002-08-06 added median keyword ; 2002-12-05 Converted to use fitsloader. ;- ;########################################################################### ; ; LICENSE ; ; This software is OSI Certified Open Source Software. ; OSI Certified is a certification mark of the Open Source Initiative. ; ; Copyright © 2001-2003 by Marshall Perrin ; ; This software is provided "as-is", without any express or ; implied warranty. In no event will the authors be held liable ; for any damages arising from the use of this software. ; ; Permission is granted to anyone to use this software for any ; purpose, including commercial applications, and to alter it and ; redistribute it freely, subject to the following restrictions: ; ; 1. The origin of this software must not be misrepresented; you must ; not claim you wrote the original software. If you use this software ; in a product, an acknowledgment in the product documentation ; would be appreciated, but is not required. ; ; 2. Altered source versions must be plainly marked as such, and must ; not be misrepresented as being the original software. ; ; 3. This notice may not be removed or altered from any source distribution. ; ; For more information on Open Source Software, visit the Open Source ; web site: http://www.opensource.org. ; ;########################################################################### ;FIXME add error checking on inputs here. if keyword_set(fdir) then datadir=fdir ; back compatibility if ~(keyword_set(datadir)) then datadir = './data/' if ~(keyword_set(outdir)) then outdir="./" file_mkdir,outdir ; this creates it if necessary checkdir,datadir ; flatfield ;;if not(keyword_set(flatf)) then $ ;; flatf = '../calib/kstwiflat451.fits' ; bad pixel mask if not(keyword_set(badpixf)) then badpixf ='' if not(keyword_set(outdir)) then cd,current=outdir ; this stores the current directory in outdir checkdir,outdir if not(keyword_set(subshift)) then subshift="S" ; default is cubic interpolation. ; FFT would work better for nyquist data but ; IRCAL is often undersampled... if not(keyword_set(submethod)) then submethod="F" ; Fourier cross-correlation fit by Gaussian if not(keyword_set(pclip)) then pclip = 0.8 ; discard the worst 20% by FWHM if ~(keyword_set(instrument)) then instrument="IRCAL" ;-------------------------------------------------- ;SDIR = '~/idl/ircal/companions/' ;SNAME = 'redircalsub.pro' getmyname,SNAME,SDIR loadct, 15 ; history preface hstr = strarr(4) hstr[0] = '<>' hstr[1] = ' istart = '+strc(istart,/join) hstr[2] = ' iend = '+strc(iend,/join) hstr[3] = ' datadir = '+strc(datadir) ; get images ; FIXME add option here to let user manually edit image list ;listims, istart, iend, out = outdir+objn+'-in.lst', $ ; datadir = datadir, fname = 'ircal', /ircal, /silent ;loadimf, outdir+objn+'-in.lst', imgs, head0, histstart = hstr, /ircal message,/info,"Now using FITSLOADER to load, bias, and flat." fitsloader,transpose([[istart],[iend]]),instr=instrument,datadir=datadir,caldir=caldir,$ imgs0,hdrs,/dobias,/doflat,listfile=listfile,badpixmask=badpix,/domask,$ flatfilename=flatf,badpixfilename=badpixf,scaletime=scaletime,$ saturation=saturation,ircam=objn imgs=temporary(imgs0) ; don't modify the originals... head0=hdrs[*,0] sxaddhist, ' REDIRCALSUB.PRO now running', head0 sxaddhist, ' Input dir: '+datadir, head0 sxaddhist, ' Input files: '+strc(istart,/join)+" - "+strc(iend,/join) , head0 sz=size(imgs) nimgs = sz[3] if ~(keyword_set(auto)) then begin print, "Examine the imgs variable now" stop endif ; output file contains exptime, rounded to the second if over 1s. time = strc(fix(sxpar(head0,"COADDTIM")*1.0/1000)) if (time eq '0') then time = strc(total(sxpararr(hdrs,"COADDTIM")*1.0/1000)) objn = objn + "."+strc(sxpar(head0,"FILTER1")) if keyword_set(verbosename) then begin objn +="-"+ time+"s"+"-"+strc(istart[0]) endif if (n_elements(badpix) eq 1) then begin $ print, '** bad pixel image does not exist !! **^G' & $ print,n_elements(badpix),(n_elements(badpix) eq 1) & $ endif else $ if n_elements(badpixf) ne 0 then sxaddhist, ' loaded bad pixel mask = "'+badpixf+'"', head0 exptimes= sxpararr(hdrs,"EXPTIME") if keyword_set(medsky) then begin ; subtract median sky frame print, '==============================' print, '=== Subtracting Median Sky ===' print, '==============================' medarr2,imgs,medsky,exptimes=exptimes display2,medsky,title="Median sky frame" wshow print, "You can examine the medsky variable now; when you press .cont the subtraction will happen" if not(keyword_set(auto)) then stop skysub,imgs,medsky,imgsout,/save,exptimes=exptimes print, " Now you can look at the imgsout variable if you want" if not(keyword_set(auto)) then stop imgs = temporary(imgsout) sxaddhist, 'SKYSUB : subtracted median sky frame', head0 endif if keyword_set(skyfile) then begin print," === Subtracting sky from "+skyfile+" ===" medsky = readfits(skyfile) if n_elements(medsky) lt 1000 then message,"Unable to read sky from "+skyfile+"!!" sxaddhist, 'SKYSUB : subtracted sky file:', head0 sxaddhist, 'SKYSUB : '+skyfile, head0 skysub,imgs,medsky,imgsout,/save,exptimes=exptimes imgs=temporary(imgsout) endif if keyword_set(submedians) then begin print, "=== Setting image medians to zero ====" for i = 0,nimgs-1 do imgs[*,*,i] -= median(imgs[*,*,i]) sxaddhist, 'SKYSUB : Set image medians to zero.', head0 endif if keyword_set(fixpix) then begin ircal_fixpix,imgs,badpix,/maskcorners,/stars sxaddhist, ' : ran ircal_fixpix on all images', head0 endif if keyword_set(satmask) then begin ircal_satmask,imgs,hdrs,after=4,instrument=instrument,satmask=satmask,/maskonly sxaddhist, ' : Masked out all saturated pixels; ',head0 sxaddhist, ' assuming pixels once saturated are bad for 4 frames',head0 endif if keyword_set(nan) then begin sxaddhist, ' : Set all NaN pixels to image medians. ',head0 message,/info,"Setting NaN pixels to image medians" message,/info,"This is just to make the cross-correlation faster" mm = meds(imgs) wnan = where(finite(imgs[*,*,0],/NaN) eq 1) sz = size(imgs) for i = 0,nimgs-1 do begin imgs[wn+i*sz[1]*sz[2]] = mm[i] endfor endif if keyword_set(dewarp) then begin platescale=0.040 ; new platescale to use in final images ;imgs0 = imgs imgs = ircal_dewarp(imgs,platescale=platescale,cubic=-0.5,/crop) sxaddhist,"IRCAL_DEWARP: regridded to "+strc(platescale)+"''/pixel",head0 ;badmask = ircal_dewarp(badmask,platescale=platescale,/byte,/nofluxscale,/crop) endif if keyword_set(shiftfile) then begin sxaddhist, ' Reading image shifts from file '+shiftfile,head0 message,/info,"Reading image shifts from file "+shiftfile readcol,shiftfile,px,py if n_elements(px) ne sz[3] or n_elements(py) ne sz[3] then $ message,"Wrong number of points in "+shiftfile endif else begin ; register all images against the first one sxaddhist, ' Finding image shifts automatically ',head0 ;stop subreg,imgs[*,*,0],imgs,shifts,method=submethod,headers=hdrs ; find the shifts subreg_shiftstopeaks,imgs,shifts,px,py sxaddhist, ' using SUBREG.PRO method='+submethod,head0 if (max(abs(px)) gt 1000) or (max(abs(py)) gt 1000) then begin message,/info,"Bogus coordinates found: Check for absurdly large position offsets!" print,"Try: print,px,py " print, "And try re-registering using a different method." stop endif endelse writecol,"./ircal_shifts.txt",px,py fwhm = fltarr(nimgs) if pclip lt 1.0 then begin ; start FWHM cut sxaddhist, ' Using FWHM cut; best '+strc(fix(pclip*100))+'% of images',head0 ; loop through images to adjust centers & get FWHMs print, form = '($,A," ")', 'FWHM computation: ' for i = 0, nimgs-1 do begin if (i mod 40) eq 0 then $ print, form = '($,A," ")', strc(i) j=i; no longer dividing by nstack ; compute FWHM radplotf, imgs[*, *, i], px[i], py[i], radout, ff, sky = 1e-20, /silent if ff eq 999 then message,"Unable to obtain FWHM!" fwhm[i] = ff print," image["+strc(i)+"] FWHM = "+strc(ff) endfor print ; plot results win, 0 w = where(fwhm ne 999) ;stop plot, fwhm, ps = 6, yr = [0, max(fwhm(w))], $ xtit = 'image number', ytit = 'FWHM (pix)' ; draw line for best subset as a guide ss = sort(fwhm) cc = (fwhm(ss))(round(nimgs*pclip)-1) oplot, !x.crange, [cc, cc] print print, 'keeping '+strc(pclip*100)+'% of the images' print, 'FWHM threshold = ', strc(cc) ; if (getyn("Use this threshhold?",1,noquery=auto) eq 0) then stop ; so you can mess with the settnigs if you want hh = fwhm < (max(fwhm(w)) < 15) if (n_elements(histogram(hh, bin = 0.5)) gt 1) then begin win, 1 plothist, hh, bin = 0.5 oplot, [cc, cc], !y.crange, line = 5 endif else print, 'not enough bins for fwhm histogram' ; select cut if ~(keyword_set(auto)) then atv,imgs,/bl fwhmcut = getnum('enter FWHM threshold:', cc,noquery=auto) goodlist = fwhm le fwhmcut wgood = where(fwhm le fwhmcut) ngood = total(goodlist) sxaddhist, ' filtered by FWHM, thrshold = '+strc(fwhmcut)+' pix', head0 endif else begin ; if pclip=1.0, then use all images. sxaddhist, ' NOT using FWHM cut; all images kept.',head0 goodlist = fltarr(nimgs)+1 wgood = findgen(nimgs) ngood=nimgs fwhmcut=" No cut made" endelse ;stop ;subpixel mosaic mosf, imgs, px, py, mos, mosexp, $ good = goodlist, badpix = badpix, masks=satmask, /setsky,$ subpixel=subpixel,interp_type=subshift,median=median,$ exptimes=exptimes,pieces=pieces if ~(keyword_set(auto)) then atv,[[[mos]],[[pieces]]],/bl if ~(keyword_set(auto)) then stop findmaxstar,mos,ix,iy mrecenter,mos,ix,iy,ixc,iyc,/silent,/nodisp radplotf, mos, ixc,iyc, radout, f1, sky = 1e-20, /silent whereismax, mos, x, y, maxv1, /silent print,"new mosaic: center at ["+strc(ixc)+", "+strc(iyc)+"]" print, 'new mosaic: FWHM = ', strc(f1), ', peak = ', strc(maxv1) sxaddhist, ' number of good images = '+strc(ngood), head0 sxaddhist, ' good images: '+strc(wgood,/join), head0 sxaddhist, ' new mosaic: FWHM = '+ strc(f1)+', peak = '+ strc(maxv1),head0 if not(getyn('satisfied?',noquery=auto)) then stop ; AO radial normalize findmaxstar,mos,xstar,ystar aoradnorm, mos, aonorm, xcen = xstar, ycen = ystar,auto=auto ircal_deghost,mos,mos_deghost,header=head0 aoradnorm, mos_deghost, aonorm_deghost, xcen = xstar, ycen = ystar,auto=auto ;;stop win, 0 & display2, mos, tit = 'mosaic' win, 1 & display2, imsub(mos, median(mos, 15)) < 20, tit = 'unsharp masked' win, 2 & display2, aonorm > (-1) < 5, /tv, tit = 'masked & filtered' win, 3 & display2, mos_deghost, tit = 'mosaic, deghosted' win, 4 & display2, imsub(mos_deghost, median(mos_deghost, 15)) < 20, tit = 'unsharp masked' win, 5 & display2, aonorm_deghost > (-1) < 5, /tv, tit = 'masked & filtere' print,"objn is "+objn if ~(keyword_set(auto)) then atv,mos,/bl if ~(keyword_set(auto)) then stop ; write files if getyn('save files?',noquery=auto) then begin if getyn('outdir="'+outdir+'": '+ $ 'do you want to change this?^G', 0,noquery=auto) then begin read, 'enter outdir: ', outdir endif file_mkdir,outdir checkdir,outdir ; write list of new shifts, FWHM, and good images outfile = outdir+objn+'_shifts.good.dat' if filecheck(outfile) then begin openw, unit0, outfile, /get_lun printf, unit0, '# RED.IDL: ', systime(), userid() printf, unit0, '# cut used for FWHM = ', strc(fwhmcut) printf, unit0, '#' printf, unit0, '# xpos ypos FWHM(pix) 1=good' printf, unit0, '' free_lun, unit0 printarray, [[px], [py], [fwhm], [goodlist]], $ out = outfile, /silent sxaddhist," Mosaic shifts and FWHMS written to "+outfile,head0 endif outfile = outdir+objn+'_mos_orig.fits' if filecheck(outfile) then begin writefits, outfile, mos, head0 if outfile ne outdir+objn+'_mos_orig.fits' then begin print,"=== Need to change objn variable for next two files to get right name" stop endif writefits, outdir+objn+'_expmap.fits', mosexp, head0 writefits, outdir+objn+'_mosaic.fits', mos_deghost, head0 writefits, outdir+objn+'_aonorm.fits', aonorm, head0 ; write parm file outfile = outdir+objn+'.parm' openw, unit0, outfile, /get_lun printf, unit0, '# REDIRCAL.IDL: ', systime(), userid() printf, unit0, '' printf, unit0, 'objn ', objn printf, unit0, 'istart ', strc(istart) printf, unit0, 'iend ', strc(iend) ;printf, unit0, 'nstack ', strc(nstack) printf, unit0, 'pclip ', strc(pclip) printf, unit0, 'datadir ', datadir ;printf, unit0, 'badpixf ', badpixf ;printf, unit0, 'flatf ', flatf printf, unit0, 'outdir ', outdir free_lun, unit0 endif endif ;;- ;;- ; save copy of reduction script ;;- ;if getyn('save reduction script in current directory?',noquery=auto) then begin ;;- if 0 then begin ;;- outfile = outdir+SNAME+'.'+objn ;;- print, ' current script dir: "', SDIR, '"' ;;- print, ' current IDL script: "', SNAME, '"' ;;- print, ' output file name: "', outfile, '"' ;;- if not(getyn(' is this correct?',noquery=auto)) then stop ;;- if filecheck(outfile,nocheck=auto) then begin ;;- spawn, 'cp '+SDIR+' '+outfile ;;- openu,lun,outfile,/append,/get_lun ;;- printf,lun,";--------------------------\n" ;;- printf,lun,"; Parameters passed to procedure:" ;;- printf,lun,"; objn= "+strc(objn) ;;- printf,lun,"; istart= "+strc(istart) ;;- printf,lun,"; iend= "+strc(iend) ;;- ;printf,lun,"; nstack= "+strc(nstack) ;;- printf,lun,"; pclip= "+strc(pclip) ;;- printf,lun,"; datadir= "+strc(datadir) ;;- if keyword_set(flatf) then printf,lun,"; flatf= "+strc(flatf) ;;- if keyword_set(badpixf) then printf,lun,"; badpixf= "+strc(badpixf) ;;- printf,lun,"; outdir= "+strc(outdir) ;;- ;;- printf,lun,"; " ;;- printf,lun,"; Script run by "+userid()+", "+systime() ;;- free_lun,lun ;;- ;;- endif ;;- ;;- print, 'wrote "', outfile, '" to current directory' ;;- endif ;;- end