pro mk_desat, infil, outdir, prefix, interactive=interactive, run_time=run_time, ist=ist ; ;+ ;NAME: ; mk_desat ;PURPOSE: ; Given a set of files, read the roadmaps and find where there ; is a long and short exposure. Make a composite image and remove ; the saturated pixels. Write images to an output file. ;INPUT: ; infil - an array of input file names ; outdir - the output directory where the file should be written ;OPTIONAL KEYWORD INPUT: ; interactive- If set, display the results to the screen as the ; processing occurs. ;HISTORY: ; draft version: 28 January 1992 (KTS) ; 20-Mar-92 (MDM) - Took Keith Strong's program "db_desat.pro" ; and made the procedure "mk_desat" ; 7-Apr-92 (KTS) - Put IFLAG in to avoid error if first image ; is rejected ; 14-Apr-92 (KTS) - 'PREFIX' parameter hardwired for 'sfd' ; - Added automatic week # feature ; - Better generated dark frame algorithm ; - Avoid night images algorithm ;Ver 2.0 16-Apr-92 (KTS) - Added dark frame library read ; 16-Apr-92 (MDM) - Changed file name definition to use the ; "progverno" variable to determine the extension ; - Put the prefix parameter back in (so it would work ; with old existing programs) ; 20-Apr-92 (MDM) - Added RUN_TIME parameter ; 21-Apr-92 (KTS) - Included LWA/MDM desat technique ; - Included MDM DC routine ;Ver 2.1 22-Apr-92 (MDM) - Added "roadmap" to the RD_SDA calls (speeds ; things up since it will not re-read the roadmap ; for each call) ; - Adjusted the format of the code so MDM can ; read it more easily ; - Fixed DC subtraction to use the whole image, ; not just one pixel ; - Put the call to GET_DC_IMAGE for the short ; exposure images outside the loop. ; - Moved the check for apropriate short exposure ; to before the RD_SDA (use the roadmap info) ; - Replaced RD_SDA with RD_XDA ; 23-Apr-92 (MDM) - Added "ist" input option to skip to image ; "ist" in the loop (for recovery from bombout) ;- ; if (n_elements(prefix) eq 0) then prefix='sfd' qdebug = 0 start_time = systime(1) progverno=2.1 ;program Version number iflag=0 ;Flag to prevent error if first ;image is rejected at late stage rd_roadmap, infil, roadmap, filidx=filidx if (keyword_set(interactive)) then window,0,xs=512,ys=512,retain=2 ; Find number of elements and set up any required arrays for longer ; exposures (omits 15 sec at present) ssla=where(roadmap.percentd eq 255 and gt_filtb(roadmap) ge 2 and $ gt_filtb(roadmap) le 3 and $ gt_res(roadmap) ge 1 and gt_expmode(roadmap) eq 0 and $ roadmap.explevmode ge 10 and roadmap.explevmode le 16) get_dc_image,roadmap(ssla),dci1,dark1,imap ; get dark frames fro long exposure dark1=sxt_decomp(dark1) sssa=where(roadmap.percentd eq 255 and gt_filtb(roadmap) ge 2 and $ ;MDM start adding 22-Apr-92 gt_filtb(roadmap) le 3 and $ gt_res(roadmap) ge 1 and gt_expmode(roadmap) eq 0 and $ roadmap.explevmode ge 2 and roadmap.explevmode le 7) get_dc_image,roadmap(sssa),dci2,dark2,imap2 ; get dark frames for short exposures dark2=sxt_decomp(dark2) imap2_arr = intarr(n_elements(roadmap)) imap2_arr(sssa) = imap2 ;MDM end adding 22-Apr-92 dla=mk_dset_str(filidx, ssla) tla=int2secarr(roadmap(ssla),roadmap(0)) nla=n_elements(tla) ; Read files sequentially and find nearest short exposure if (keyword_set(ist)) then iflag = 1 ;file already started if (n_elements(ist) eq 0) then ist = 0 for i = ist,nla-1 do begin time1_all = systime(1) rr = roadmap(ssla(i)) print, i+1, ' of ', nla,' images ', fmt_tim(rr), gt_res(rr,space=3) ; Print which frame we are on ss=where(roadmap.percentd eq 255 and $ gt_filtb(roadmap) eq gt_filtb(rr) and $ gt_res(roadmap) eq gt_res(rr) and $ gt_expmode(roadmap) eq 0 and $ roadmap.explevmode ge 2 and $ roadmap.explevmode le 7) if (ss(0) lt 0) then goto, lab1 ;no short exp - loop ;;rd_sda,infil(dla(i).ifil),dla(i).dset,index1,data1 ;read next long exp. ; search for short exposure to match it ;;ss=where(roadmap.percentd eq 255 and $ ;; gt_filtb(roadmap) eq gt_filtb(index1) and $ ;; gt_res(roadmap) eq gt_res(index1) and $ ;; gt_expmode(roadmap) eq 0 and $ ;; roadmap.explevmode ge 2 and $ ;; roadmap.explevmode le 7) ;; if (ss(0) lt 0) then goto, lab1 ;no short exp - loop ds=mk_dset_str(filidx,ss) ; Find the closest time of a short exposure tss = int2secarr(roadmap(ss),roadmap(0)) ;;dt = tss - tla(i) ;;ss = where(dt lt 0) ;;if (ss(0) ge 0) then dt(ss) = -dt(ss) dt = abs(tss - tla(i)) ;MDM replaced above 22-Apr-92 - need "ss" preserved for below, and replaced code is "cleaner" ctm = min(dt) if (ctm gt 600.) then goto, lab1 ;too long - loop id=where(dt eq ctm) ;id is in reference to ss subset ;;rd_sda,infil(ds(id(0)).ifil),ds(id(0)).dset,index2,data2 time1 = systime(1) rd_xda, infil, ds(id(0)), index2, data2, roadmap time2 = systime(1) if (qdebug) then print, 'RD_XDA for short exposure took ', (time2-time1), ' seconds' time1 = systime(1) rd_xda, infil, dla(i), index1, data1, roadmap ;read next long exp. time2 = systime(1) if (qdebug) then print, 'RD_XDA for long exposure took ', (time2-time1), ' seconds' ;;get_dc_image,index2,dci2,back2 ; Find a dark frame ;MDM removed 22-Apr-92 ;;back2d=sxt_decomp(back2) time1 = systime(2) datal=sxt_decomp(data1) ; decompress short and long exposures datas=sxt_decomp(data2) time2 = systime(2) if (qdebug) then print, 'Decompression for two images took ', (time2-time1), ' seconds' ; Short exposure taken at night? - if so then loop ; if (max(datas) lt 5000.) then goto, lab1 ; if (max(datal) lt 5000.) then goto, lab1 ; Subtract background and normalize exposures ;datal=datal-dark1(imap(i)) ;this was just subtracting a scalar value!! - MDM removed 22-Apr-92 ;datas=datas-back2d time1 = systime(2) datal=datal-dark1(*,*,imap(i)) datas=datas-dark2(*,*,imap2_arr(ss(id(0)))) ;MDM added 22-Apr-92 time2 = systime(2) if (qdebug) then print, 'Dark image subtraction took ', (time2-time1), ' seconds' if (qdebug) then begin temp1 = dci1(imap(i)) temp2 = dci2(imap2_arr(ss(id(0)))) print, 'Long: ' + fmt_tim(index1) + gt_res(index1,sp=2) + gt_dpe(index1,/conv,/str) + $ ' Dark: ' + fmt_tim(temp1) + gt_res(temp1,sp=2) + gt_dpe(temp1,/conv,/str) print, 'Short: ' + fmt_tim(index2) + gt_res(index2,sp=2) + gt_dpe(index2,/conv,/str) + $ ' Dark: ' + fmt_tim(temp2) + gt_res(temp2,sp=2) + gt_dpe(temp2,/conv,/str) end time1 = systime(2) datal=datal*(1000./gt_expdur(index1)) datas=datas*(1000./gt_expdur(index2)) datal=datal>1 datas=datas>1 time2 = systime(2) if (qdebug) then print, 'Image normalization took ', (time2-time1), ' seconds' time1 = systime(2) m1 = where(datas ge 90) ;find good pixels in short s2=where(data1 ge 253) ;find sat. pixels if (gt_res(index1) eq 1) then begin xx= s2 mod 512 yy = s2/512 s2up = ((yy+1)<511)*512L + xx s2dn = ((yy-1)>0 )*512L + xx endif else begin xx=s2 mod 256 yy = s2/256 s2up = ((yy+1)<255)*256L + xx s2dn = ((yy-1)>0)*256L + xx endelse ; Put the longest exposure in the working array ; and replace saturated pixels with smoothed med image ;;cimg=datal ;MDM removed and renamed all references to "cimg" to "datal" ;save memory/computing time if (s2(0) ne -1) then begin ;MDM added 22-Apr-92 (some images have no saturated pixels) datal(s2) = datas(s2) datal(s2up) = datas(s2up) datal(s2dn) = datas(s2dn) timg=smooth(datal,5) datal(s2) = timg(s2) datal(s2up) = timg(s2up) datal(s2dn) = timg(s2dn) end if (m1(0) ne -1) then datal(m1)=datas(m1) ; Replace "m1" pixels in datal with unsmoothed med exposure values ;;datal=cimg ;Keith, why do this copy? why not use "cimg" below? (MDM 22-Apr-92) ; normalize for different resolutions ;;factor=1. ;;if (gt_res(index1) eq 2) then factor=factor*0.25 ;;datal=factor*datal if (gt_res(index1) eq 2) then datal=datal*0.25 ;MDM added - avoid extra CPU math when possible ; compress the images to byte arrays print,'Maximum Count Rate: ',max(datal) datal=alog10(datal>1.) datal=bytscl(datal,max=6.,min=0.,top=255) print,'Maximum Count Rate: ',max(datal) time2 = systime(2) if (qdebug) then print, 'Desaturation/image combination/image scaling took ', (time2-time1), ' seconds' if (keyword_set(interactive)) then begin if (gt_res(index1) eq 2) then begin tv,rebin(datal,512,512,/sample) endif else begin tv,datal endelse text= string(i) + ' ' + fmt_tim(index1) + ' ' + $ gt_filtb(index1,/string) + ' ' + gt_res(index1,/string) xyouts,10,10,text,/device end time1 = systime(2) if ((iflag eq 0) or (i eq ist)) then begin iflag=1 ;reset flag anum=strlen(infil(0)) ;length of file id fileid = getobsid(startt=strmid(infil(0),anum-11,6)) ;create Week ID fileid=fileid + '.' + string(fix(progverno), format='(i2.2)') filename = concat_dir(outdir, prefix+fileid) end if (iflag eq 0) then begin sav_sda,filename,index1,datal end else begin sav_sda,filename,index1,datal,/append end time2 = systime(2) if (qdebug) then print, 'SAV_SDA took ', (time2-time1), ' seconds' time2_all = systime(1) if (qdebug) then print, '*********************** One image process time = ', (time2_all-time1_all), ' seconds lab1: dumy=0 endfor end_time = systime(1) run_time = (end_time-start_time)/60. end