;****************************************************************************** ;Given a list of BCD images, return a coadded image. ; ;INPUTS: ; ;imagelistname = required list of BCD images to coadd, one image name per ;row. There must be more than one image in the list, or a warning ;will be issued and the program will not create any outputs. ; ;verbose = optional keyword. If set, information about the processing ;will be printed to the screen. ; ;fatalbits = optional keyword allowing the user to set the bits ;considered to be fatal. These will be set to NaN in the BCD and ;uncertainty files before coaddition. The input keyword should be an ;array containing the fatal bits. The default is ;fatalbits = [8, 12, 13, 14]. If a user wishes to have no fatal ;bits, they should set fatalbits = [999]. ; ;OUTPUTS: ; ;coadded image - *coa2d.fits ;uncertainty image - *c2unc.fits ;mask file - *c2msk.fits ; ;EXAMPLE: ; ;coad, 'bcdlist.txt', /verbose, fatalbits=[8,12,13,14] ; ;****************************************************************************** pro coad, imagelistname, verbose=verbose, fatalbits=fatalbits if keyword_set(fatalbits) then fatal_bits = fatalbits else fatal_bits = [8, 12, 13, 14] if keyword_set(verbose) then begin print, '' print, 'Bits in bmask files set as fatal: ', fatal_bits print, '' endif ;Read in the list of BCD images to combine. readcol, imagelistname, imagelist, format = '(a)', /silent num = n_elements(imagelist) if keyword_set(verbose) then begin print, '' print, 'Number of images to combine: ', num print, '' print, imagelist print, '' endif ;Print out a warning if there is only one image to combine, and ;don't do anything else. if num gt 1 then begin ;Figure out the names of the corresponding uncertainty files. unc_imagelist = strarr(num) for i = 0, num - 1 do begin unc_imagelist[i] = repstr(imagelist[i], 'bcd.fits', 'func.fits') endfor if keyword_set(verbose) then begin print, 'Uncertainty files: ' print, '' print, unc_imagelist print, '' endif ;Figure out the names of the corresponding bmask files. mask_imagelist = strarr(num) for i = 0, num - 1 do begin mask_imagelist[i] = repstr(imagelist[i], 'bcd.fits', 'bmask.fits') endfor if keyword_set(verbose) then begin print, 'Mask files: ' print, '' print, mask_imagelist print, '' endif ;Figure out the output file name based on the first input BCD name. nameparts = strsplit(imagelist[0], '_', /extract) outfilename = nameparts[0]+'_'+nameparts[1]+'_'+nameparts[2]+'_'+nameparts[3]+'_'+nameparts[5]+'_coa2d.fits' channel = repstr(nameparts[1], 'S', '') ;How big is each image? first_image = readfits(imagelist[0], hdr, /silent) s = size(first_image) xsize = s[1] ysize = s[2] ;Initialize the data cube. cube = fltarr(xsize, ysize, num) unc_cube = fltarr(xsize, ysize, num) ;Fill the data cube. for i = 0, n_elements(imagelist) - 1 do begin ;Read in the BCD, uncertainty, and mask images. image = readfits(imagelist[i], hdr, /silent) unc_image = readfits(unc_imagelist[i], unc_hdr, /silent) mask_image = readfits(mask_imagelist[i], mask_hdr, /silent) ;Set all pixels with fatal bits to NaN in bcd and uncertainty files. for ix = 0, xsize - 1 do begin for iy = 0, ysize - 1 do begin for b = 0, n_elements(fatal_bits) - 1 do begin flag = bitget(mask_image[ix,iy], fatal_bits[b]) if flag eq 1 then begin image[ix,iy] = !Values.F_NAN unc_image[ix,iy] = !Values.F_NAN endif endfor endfor endfor cube[0,0,i] = image unc_cube[0,0,i] = unc_image endfor ;Compute the mask file. mask_image = readfits(mask_imagelist[0], mask_hdr, /silent) for i = 1, num - 1 do begin nextmask = readfits(mask_imagelist[i], temp_hdr, /silent) mask_image = mask_image or nextmask endfor ;Compute the trimmed mean of the data cube and its uncertainty. output_image = dblarr(xsize, ysize) unc_image = dblarr(xsize, ysize) for i = 0, xsize - 1 do begin for j = 0, ysize - 1 do begin ;mask file. bit7_flag = bitget(mask_image[i,j], 7) if bit7_flag eq 1 then mask_image[i,j] = 2L^7 ;trimmed mean data = reform(cube[i,j,*]) unc_data = reform(unc_cube[i,j,*]) finite_flag = finite(data) finite_index = where(finite_flag eq 1, finite_count) if finite_count eq 0 then begin output_image[i,j] = !Values.F_NAN mask_image[i,j] = mask_image[i,j] + 2L^13 endif if finite_count eq 1 then output_image[i,j] = data[finite_index] if finite_count ge 2 then begin meanclip, data[finite_index], trimmed_mean output_image[i, j] = trimmed_mean endif ;uncertainty unc_finite_flag = finite(unc_data) unc_finite_index = where(unc_finite_flag eq 1, unc_finite_count) if unc_finite_count eq 0 then unc_image[i,j] = !Values.F_NAN if unc_finite_count gt 0 then begin unc_image[i, j] = sqrt(total((unc_data[unc_finite_index])^2.)) / double(unc_finite_count) endif endfor endfor if keyword_set(verbose) then begin print, 'Computing trimmed mean of input BCD files.' print, '' print, 'Computing output uncertainty as sum of squares of input BCD uncertainty files, normalized by number of valid images for a given pixel.' print, '' endif ;Write output files. writefits, outfilename, float(output_image), hdr writefits, repstr(outfilename, 'coa2d.fits', 'c2unc.fits'), float(unc_image), unc_hdr writefits, repstr(outfilename, 'coa2d.fits', 'c2msk.fits'), mask_image, mask_hdr if keyword_set(verbose) then begin print, 'Writing output coad file: ', outfilename print, '' print, 'Writing output uncertainty file: ', repstr(outfilename, 'coa2d.fits', 'c2unc.fits') print, '' print, 'Writing output mask file: ', repstr(outfilename, 'coa2d.fits', 'c2msk.fits') endif endif else begin print, '' print, 'WARNING: Only one file name in '+imagelistname+'. No coad created.' print, '' endelse end