;+ ; NAME: ; cube2biasdark ; ; PURPOSE: (one line) ; make zero-flux (dark+bias) averaged image from image cube ; ; DESCRIPTION: ; ; ; CATEGORY: ; IMAGES ; ; CALLING SEQUENCE: ; biasdark=cube2dark(cube, head, err, bad, d = d, $ ; silent=silent, median = median, $ ; imbad = imbad, bias=bias, minbias=minbias, $ ; readn_dn = 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_zero1.pro ; Written 14 Jun 2007 Leslie Ann Young SwRI ; Modified 13 Sep 2007 LAY ; 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 cube2dark ; 2010-08-30 --- convert from make_spe_zero1 ;- pro cube2dark, cube, avg, err, bad, header=header, readn_dn = readn_dn, $ silent=silent, median = median, $ imbad = imbad, in_readn_dn = in_readn_dn 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 ;-------------------------------- ; Set readnoise (rn). ; If the readnoise is not specified, estimate it from the ; difference in the first two images ;-------------------------------- if keyword_set(in_readn_dn) then begin rn = in_readn_dn 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 rn = stddev endif else begin d0 = float(d[*,*,1]) - float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal rn = stddev / sqrt(2.) endelse endelse if not silent then begin print, 'READNOISE (1st pass) = ', rn, 'DN' end ;-------------------------------- ; Get the zero/bias frame ; ;-------------------------------- if keyword_set(median) then method = 'median' else method = 'avgclipsim' case nimages of ;-------------------------------- ; 1 image: ; avg = that image ; err = passed read noise, or estimated read noise 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(rn, nx, ny) ; error in the zero smavg = median(avg, 10) bad = bytarr(nx,ny) indxbad = where( abs(avg - smavg)/err gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end ;-------------------------------- ; 2 images: ; avg = mean of the two images ; err = passed read noise, or estimated read noise 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(rn, nx, ny)/2. ; error in the zero bad = bytarr(nx,ny) indxbad = where (abs(d[*,*,0] - d[*,*,1])/(sqrt(2.)*rn) gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end 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 print, 'READNOISE = ', avgerr, 'DN' ; 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, avgclip algorithm sigma = rn badmat = bytarr(nx,ny,nimages) 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] resid = (new-avg)/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 = replicate(rn, nx, ny) 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 if not keyword_set(silent) then begin print, 'READNOISE = ', avgerr, 'DN' endif end ; avgclipsim endcase end end readn_dn = avgerr end