pro irac_subcube_collapse, ifile, ofile, sfile, INPUT_MASK=input_mask, $ NSIGMA=nsigma, COUNT_IMAGE=count_image, $ OUTPUT_MASK=output_mask ;+ ; NAME: ; IRAC_SUBCUBE_COLLAPSE ; ; PURPOSE: ; Procedure takes an input IRAC subarray cube BCD file and generates a ; corresponding two dimensional analog that can be used to generate mosaics ; and perform source extraction using the SSC Mopex software. The ; procedure also generates a corresponding uncertainty image using the ; variance of the pixels at each array location that are used to determine ; the final pixel value. ; ; INPUT: ; ifile: filename of subarray BCD to convert ; (e.g. SPITZER_I1_13289216_0008_0000_4_bcd.fits ; ; OUTPUT: ; ofile: filename of 2d dimensional version of asubarray BCD ; sfile: filename of corresponding 2d uncertainty file ; ; OPTIONAL INPUT KEYWORDS: ; INPUT_MASK: input dmask, if used, then output mask is OR of mask pixels ; that go into collapsed pixel stacks ; NSIGMA: Number of sigma rms used in outlier rejection for trimmed means. ; If not set, then a default of 3 is used. ; ; OPTIONAL OUTPUT KEYWORD: ; COUNT_IMAGE: if set, then write image containing the number of subarray ; samples used at each position of the output BCD ; OUTPUT_MASK: if set (and if INPUT_MASK) is set, then an output mask which ; is the OR of all the good mask pixels in the stack is made ; ; ; METHOD: ; The input BCD cube (nx by ny by nz samples) is read in. nx by ny pixel ; output data and uncertainty images are created. For each i, j pixel ; position, the output image value is the outlier trimmed mean of the ; input pixel stack I[i, j, *]. N sigma iterative rejection is used ; with the outliers being recalculated until the number of rejected pixels ; in the stack remains constant. N is set to 3 by default, but is ; controllable by the NSIGMA keyword. The corresponding uncertainty value of ; for pixel i, j is the variance of the pixels used to calculate the mean. ; An optional output image which is the number of samples from the stack ; used for each output pixel is produced if the COUNT_IMAGE keyword is ; set. If INPUT_MASK and OUTPUT_MASK are set, then a mask image is ; generated which is the OR of all the mask values in the stack of pixels ; used in the final image. A header is geneated for the output files from ; the input cube (and mask) header. The third dimension is stripped from ; the header keywords. The two dimensional image, uncertainty and mask ; are written to the specified filenames. ; ; FUNCTIONS USED: ; READFITS, MOMENT, FILE_SEARCH, SXPAR ; ; PROCEDURES USED: ; WRITEFITS, SXADDPAR ; ; HISTORY: ; Code adopted from S. Fajardo-Acosta by S. J. Carey 06 Nov 2006 ; Corrected bug which cause routine to make only one pass through a stack ; S. J. Carey 27 Nov 2007 ; Error in uncertainty calculation, did not divide by sqrt(number of images ; in stack) S. J. Carey 28 Nov 2007 ;- ; Give execution information if (N_params() lt 2 or N_params() gt 3) then begin print, 'Syntax -- IRAC_SUBCUBE_COLLAPSE, incubefile, out_2dfile, $' print, ' out_2dsig, NSIGMA=nsigma, COUNT_IMAGE=cfile, $ ' print, ' INPUT_MASK=input_mask, OUTPUT_MASK=output_mask' print, ' ----------------------------------------------------------' print, ' incubefile is the subarray BCD to collapse' print, ' out_2dfile is the resulting 2d image' print, ' out_2dsig is the resulting 2d uncertainty image' print, ' out_2dmask is the resulting 2d mask (optional)' print, ' NSIGMA if set is outlier threshold' print, ' COUNT_IMAGE if set is 2d coverage map' print, ' INPUT_MASK if set is 3d imask or dmask' print, ' OUTPUT_MASK if set and INPUT_MASK is set, will' print, ' contained 2d version of INPUT_MASK' print, ' ----------------------------------------------------------' print, ' procedure assumes that you are collapsing the cube' print, ' about the third axis' return endif ; Set sigma threshold if (not keyword_set(NSIGMA)) then nsigma = 3.0 ; Check input file if (file_search(ifile) eq '') then begin print, 'IRAC_SUBCUBE_COLLAPSE: ' + ifile + ' not found' return endif ; Check input mask if (keyword_set(INPUT_MASK)) then begin if (file_search(input_mask) ne '') then begin do_mask = 1 mask = readfits(input_mask, hmask, /SILENT) ; Check to make sure it is a cube naxis = sxpar(hmask, 'NAXIS') sz = size(mask) if (naxis ne 3 or sz[0] ne 3) then begin print, 'IRAC_SUBCUBE_COLLAPSE: ' + ifile + ' is not a FITS cube' return endif ; Prepare output mask sxaddpar, hmask, 'NAXIS', 2 ; Remove third axis sxdelpar, hmask, 'NAXIS3' ; Copy header info to sigma header endif else begin print, 'IRAC_SUBCUBE_COLLAPSE: ' + input_mask + ' not found' do_mask = 0 endelse endif else do_mask = 0 ; Read in input cube cube = readfits(ifile, hcube, /SILENT) ; Check to make sure it is a cube naxis = sxpar(hcube, 'NAXIS') sz = size(cube) if (naxis ne 3 or sz[0] ne 3) then begin print, 'IRAC_SUBCUBE_COLLAPSE: ' + ifile + ' is not a FITS cube' return endif ; Going to assume that the fits file is properly formatted at this stage ; Copying header values to ouput header himage = hcube sxaddpar, himage, 'NAXIS', 2 ; Remove third axis sxdelpar, himage, 'NAXIS3' ; Copy header info to sigma header hsigma = himage nx = sxpar(himage, 'NAXIS1') ny = sxpar(himage, 'NAXIS2') ; Now create output images -- initialize to not a number image = fltarr(nx, ny) + !VALUES.F_NAN sigma = fltarr(nx, ny) + !VALUES.F_NAN omask = intarr(nx, ny) ; Number of pixels from stack per position icount = intarr(nx, ny) ; Now determine the filtered mean and variance for each (x, y) in the cube ; Loop over all x ...., hate to loop but for IRAC images this should be fast for i = 0, nx-1 do begin ; Loop over all y .... for j = 0, ny-1 do begin ; Initialize counters ; Allow a maximum of 10 iterations maxiter = 10 iter = 1 ; Number of pixels in stack from previous iteration ocount = 0 ; Grab a pixel stack vec = cube[i, j, *] ; Remove NaNs and find number of good pixels in stack ptr = where(vec eq vec, count) if (count gt 0) then vec = vec[ptr] ; If there is an input mask, get that pixel stack as well if (do_mask eq 1) then begin mvec = mask[i, j, *] if (count gt 0) then mvec = mvec[ptr] endif ; While the number of pixels in the stack is changing and we haven't reached ; the max number of iterations and there are still enough pixels to calculate ; a first moment, performing sigma rejection on stack while (count ne ocount and iter lt maxiter and count gt 2) do begin mom = moment(vec) vsig = sqrt(mom[1]) v0 = mom[0] ; Error in setting old counter, should be done before I reset the counter, and ; it will be changed on the next iteration through the loop ocount = count ptr = where(abs(vec - median(vec)) le nsigma * vsig, count) ; If there are outliers above a threshold, then trim stack if (count ne ocount) then begin vec = vec[ptr] if (do_mask eq 1) then mvec = mvec[ptr] iter = iter + 1 endif endwhile ; Place mean value and uncertainty in output images if (count gt 2) then begin image[i, j] = v0 ; Originally using standard deviation of stack as uncertainty, but as ; final value is the mean of the stack, the uncertainty should be the coadded ; standard deviation SJC 28 Nov 2007 ; sigma[i, j] = vsig sigma[i, j] = vsig / sqrt(count) icount[i, j] = count ; If input mask is used, then final mask is OR of all mask pixels remaining in ; stack if (do_mask eq 1) then for k = 0, count-1 do $ omask[i, j] = omask[i, j] or mvec[k] endif endfor endfor ; All finished write output files sxaddpar, himage, 'COLCUBE', 'T', ' Collapsed 3d cube' sxaddpar, hsigma, 'SIGCUBE', 'T', ' Sigma of Collapsed 3d cube' hcount = himage sxaddpar, hcount, 'BUNIT', 'Samples', ' Number of samples in output image' if (do_mask eq 1) then $ sxaddpar, hmask, 'COLMASK', 'T', ' Mask is the collapsed dmask' writefits, ofile, image, himage writefits, sfile, sigma, hsigma if (keyword_set(OUTPUT_MASK)) then writefits, output_mask, omask, hmask if (keyword_set(COUNT_IMAGE)) then writefits, count_image, icount, hcount return end