;+ ; NAME: ; cube2flat ; ; PURPOSE: (one line) ; make flat-frame averaged image from image cube ; ; DESCRIPTION: ; ; ; CATEGORY: ; IMAGES ; ; CALLING SEQUENCE: ; biasdark = cube2dark, cube, avg, err, bad, header=header, readn_dn = readn_dn, $ ; silent=silent, median = median, $ ; imbad = imbad, in_readn_dn = in_readn_dn ; INPUTS: ; cube : ; ; OPTIONAL INPUT PARAMETERS: ; -------- behavior of the routine ------------ ; silent : 1 if the routine is to perform silently. Default = 0, ; print out statements ; median : 1 if using median averaging, 0 if using avgclip algorithm ; ; -------- description of the input file ------ ; imbad : if set, an array of bad image indices (eg ; imbad=[0] to drop 1st image ; set_readn_dn are used in some modes to calculate the ; expected noise per pixel, used for flagging outliers ; ; OUTPUTS: ; avg - the average zero-flux image ; err - the error in the average zero-flux image (NOT error per ; image per pixel) ; bad - bad pixel image for avg (0 = good) ; readn_dn : read noise in dn. Default = caculate from the data ; ; OPTIONAL OUTPUTS: ; d - the input file (after deselecting with imbad) ; ; EXAMPLE: ; ; MODIFICATION HISTORY: ; HISTORY OF make_spe_flat1.pro ; Written 14 Jun 2007 Leslie Ann Young SwRI ; Modified 13 Sep 2007 LAY ; Modified 18 Sep 2007 LAY added exptime keyword ; Modified 20 Sep 2007 LAY print input file ; Modified 21 Sep 2007 LAY median only the middle 1000 frames ; Modified 22 Sep 2007 CBO fixed two minor bugs from median'ing the middle 1000 frames ; HISTORY OF cube2flat ; 2010-09-22 --- convert from make_spe_flat1 ;- pro cube2flat, cube, avg, err, bad, header=header, gain = gain, $ zero = zero, zerr=zerr, zbad=zbad, znoise=znoise, $ silent=silent, median = median, $ imbad = imbad, in_gain = in_gain if not keyword_set(silent) then silent = 0b d = cube nxnyn, d, nx, ny, nimages ;-------------------------------- ; Filter out the bad images ;-------------------------------- if n_elements(imbad) gt 0 then begin if not silent then begin print, 'DESELECTING IMAGES ', imbad end isgood = bytarr(nimages) + 1b isgood[imbad] = 0 igood = where(isgood eq 1b, ngood) if ngood eq 0 then begin print, 'make_spe_zero1(', filename,'...): all images deselected' return endif else begin d = d[*,*, igood] endelse nxnyn, d, nx, ny, nimages end ;-------------------------------- ; Subtract the zero-flux ;-------------------------------- if not keyword_set(zero) then zero = fltarr(nx, ny) if not keyword_set(zerr) then zerr = fltarr(nx, ny) if not keyword_set(zbad) then zbad = bytarr(nx, ny) if not keyword_set(znoise) then znoise = 0. d = float(d) for i = 0, nimages-1 do begin d[*,*,i] = d[*,*,i] - zero endfor ;-------------------------------- ; Set gain ; If the gain is not specified, estimate it from the ; first image ;-------------------------------- if keyword_set(in_gain) then begin g = in_gain endif else begin if nimages eq 1 then begin d0 = float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal err = stddev endif else begin d0 = float(d[*,*,1]) - float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal ; get noise err = stddev / sqrt(2.) d0 = float(d[*,*,1]) + float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal ; get level endelse g = avg/(err^2 - znoise^2) endelse if not silent then begin print, 'GAIN (1st pass) = ', g, 'e/DN' end ;-------------------------------- ; Get the flat field frame ; ;-------------------------------- if keyword_set(median) then method = 'median' else method = 'avgclipsim' means = fltarr( nimages ) for i=0,nimages-1 do begin skysclim,d[*,*,i],lowval,hival,meanval,sigma,npts=20000 means[i] = meanval d[*,*,i]=d[*,*,i]/means[i] if not silent then begin if i lt 20 then print,'Frame ',i,' scaled by ',means[i] endif endfor normfac = median(means) case nimages of ;-------------------------------- ; 1 image: ; avg = that image ; err = passed gain, or estimated gain from robomean ; bad = 1 where points are more than 5 sigma from a 10x10 median-smoothed ; version of the average ;-------------------------------- 1: begin thresh = 5.0 avg = d[*,*,0] err = replicate( sqrt(means[0]/g), nx, ny) ; error in the flat smavg = median(avg, 10) bad = bytarr(nx,ny) if keyword_set(inbadmat) then begin bad = inbadmat[*,*,0] endif indxbad = where( abs(avg - smavg)/err gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end ; nimage eq 1 ;-------------------------------- ; 2 images: ; avg = mean of the two images ; err = passed gain, or estimated gain from robomean ; bad = 1 where the difference is > 5 sigma (5 * rn * sqrt(2)) ;-------------------------------- 2: begin thresh = 5. avg = (d[*,*,0] + d[*,*,1])/2. err = replicate(sqrt((means[0]+means[1])/g), nx, ny)/2. ; error in the zero bad = bytarr(nx,ny) if keyword_set(inbadmat) then begin d0 = d[*,*,0] d1 = d[*,*,1] b0 = inbadmat[*,*,0] b1 = inbadmat[*,*,1] ; the points that are bad in frame 0 and good in frame 1 bad0 = where(b0 eq 1 and b1 eq 0, nbad0) if nbad0 gt 0 then begin avg[bad0] = d1[bad0] err[bad0] = sqrt(2.) * err[bad0] endif ; the points that are good in frame 0 and bad in frame 1 bad1 = where(b0 eq 0 and b1 eq 1, nbad1) if nbad1 gt 0 then begin avg[bad1] = d0[bad1] err[bad1] = sqrt(2.) * err[bad1] endif ; the points that are bad in both bad01 = where(b0 eq 1 and b1 eq 1, nbad1, compl=good01, ncompl=ngood01) if ngood01 eq 0 then begin print,' Serious Warning! all pixels in imge cube were flagged bad.' bad[*,*] = 1 endif else begin if nbad1 gt 0 then begin avg[bad1] = median(avg[good01]) bad[bad01] = 1 endif endelse endif indxbad = where (abs(d[*,*,0] - d[*,*,1])/(sqrt(2.)*rn) gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end ; nimage eq 2 else: begin case method of ;-------------------------------- ; >=3 images, median: ; avg = median ; err = std deviation of mean sqrt( total(resid^2)/( (n-1)*n) ) ; bad = 1 where error is 5 sigma larger than other errors ;-------------------------------- 'median': begin ; call the (slow) medarr_mwb rather ; than risk malloc error if nimages gt 1000 then begin medarr_mwb, d, avg endif else begin avg = median(d, dim=3, /even) endelse ; create the error array res = d for i = 0, nimages-1 do res[*,*,i] = d[*,*,i]-avg ; err1 is the error in one pixel (i.e., the readnoise) err1 = fltarr(nx, ny) for i=0, nx-1 do begin for j = 0, ny-1 do begin err1[i,j] = total(res[i,j,*]^2) endfor endfor err1 = sqrt( err1/(nimages-1) ) ; error in a single pixel err = err1/sqrt(nimages) ; error in the average robomean,err1,5.,0.1,avgerr,avgerrdev,stddeverr gnew = (avgerr^2 - (znoise/normfac)^2)/normfac if gnew ge 0. then g = gnew if not keyword_set(silent) then begin print, 'GAIN = ', g, 'e/DN' endif ; create the bad-pixel array bad = bytarr(nx, ny) ibad = where(err gt (avgerr + 5.0 * stddeverr), nbad) if nbad gt 0 then bad[ibad] = 1 end ; median ;-------------------------------- ; >=3 images, median: ; avg = average, clipping 3-sigma outliers ; err = std deviation of mean excluding marked outliers ; bad = 1 where all pixels are outliers ;-------------------------------- 'avgclipsim': begin thresh = 3.0 ; 2007 Sept 21: malloc errors when ; median'ing 0.5-Gb file if nimages gt 1000 then begin first = (nimages / 2 - 500) > 0 last = (first + 999) < (nimages - 1) avg = median(d[*,*, first:last], dim=3, /even) endif else begin avg = median(d, dim=3, /even) endelse ; approximate noise sigma = sqrt( ( (avg>0)*normfac/g) + znoise^2 + zerr^2 ) $ / normfac if keyword_set(inbadmat) then begin badmat = inbadmat endif else begin badmat = bytarr(nx,ny,nimages) endelse asum = fltarr(nx,ny) acnt = fltarr(nx,ny) npts = long(nx) * long(ny) npixclipped = 0 for i=0,nimages-1 do begin new = d[*,*,i] newavg = avg * median(new)/median(avg) resid = (new - newavg)/sigma z = where(abs(resid) ge thresh,count) if count ne 0 then begin badi = badmat[*,*,i] badi[z] = 1B badmat[*,*,i] = badi endif z = where(badmat[*,*,i] eq 0B,count) if count ne 0 then begin asum[z] = asum[z]+new[z] acnt[z] = acnt[z] + 1.0 endif else begin print,' Serious Warning! all pixels in image ',strn(i), $ ' were considered bad.' ;stop endelse npixclipped = npixclipped + (npts - count) endfor ; refine the average where it's possible to z = where(acnt ne 0,count) if count ne 0 then avg[z]= asum[z]/acnt[z] ; create the bad pixel array bad = (acnt le 1) ; create the array with the error in avg err1 = sigma for i = 0, nx-1 do begin for j=0, ny-1 do begin n = acnt[i,j] if n ge 2 then begin err1[i,j] = $ sqrt( total( $ (d[i,j,*] - avg[i,j])^2 * (1.-badmat[i,j,*])) / $ (n-1.) ) endif endfor endfor err = err1 / sqrt(nimages) robomean,err1,5.,0.1,avgerr,avgerrdev,stddeverr gnew = normfac/( (avgerr*normfac)^2 - znoise^2 - median(zerr)^2 ) if gnew ge 0. then g = gnew if not keyword_set(silent) then begin print, 'GAIN = ', g, 'e/DN' endif end ; avgclipsim endcase ; method end ; nimage gt 3 endcase ; nimages gain = gnew end