;create mask image according to mask file ;astrometry keyword for image input use to convert ra,dec mask ;coordinates into pixel using the image WCS ;Structure obtained from extast.pro procedure reading the header of a ;fits file ; pro deg2radec_2d,ra_cent, dec_cent, x, y, ra_out, dec_out x_t = reform(tan(x*!dpi/180.)) y_t = reform(tan(y*!dpi/180.)) x0 = ra_cent*(!dpi/180.) y0 = dec_cent*(!dpi/180.) DET = cos(y0)-y_t*sin(y0) size_x_t = size(x_t, /dimension) size_y_t = size(y_t, /dimension) ra_out = dblarr(size_x_t[0], size_x_t[1]) dec_out = dblarr(size_y_t[0], size_y_t[1]) ra_out = (x0 -atan(x_t,DET))*180./!dpi z = (cos(ra_out *!dpi/180.-x0))/DET dec_out = (atan(z*(sin(y0) +y_t*cos(y0))))*180./!dpi end function radec2xyz, ra, dec deg2rad = !dpi/180 ra *= deg2rad dec *= deg2rad cos_ra = cos(ra) cos_dec = cos(dec) sin_ra = sin(ra) sin_dec = sin(dec) n = size(reform(ra), /dimension) ;array case res = dblarr(3,n[0], n[1]) res[0,*,*] = cos_ra * cos_dec res[1,*,*] = sin_ra * cos_dec res[2,*,*] = sin_dec return, res end function get_cxyz_image, sx, sy, ra_min, dec_min, spix, astrometry = $ astrometry, radec_cent = radec_cent ;get ra, dec on sky if(keyword_set(astrometry)) then begin ;get ra, dec for image pixels x_arr = lindgen(sx) y_arr = lindgen(sy) xy2ad, x_arr, y_arr, astrometry, ra_arr, dec_arr ;build 2d image radec_img = dblarr(2,sx, sy) ;ra radec_img[0,*,*] = cmreplicate(ra_arr, [sy]) ;dec radec_img[1,*,*] = transpose(cmreplicate(dec_arr, [sx])) endif else begin ;get ra_min, dec_min on SKY ;radec_min_sky = deg2radec(radec_cent[0], radec_cent[1], ra_min, dec_min) x = ra_min + spix/2. + spix*dindgen(sx) y = dec_min + spix/2. + spix*dindgen(sy) ;build 2d image x_img = cmreplicate(x, [sy]) y_img = transpose(cmreplicate(reverse(y), [sx])) deg2radec_2d, radec_cent[0], radec_cent[1], x_img, y_img, ra_img, dec_img radec_img = dblarr(2,sx, sy) ;reverse because will be used on ;equator projected coordinates, which ;here, have ra reversed radec_img[0,*,*] = reverse(ra_img) radec_img[1,*,*] = dec_img endelse ;get cxyz cxyz = radec2xyz(reform(radec_img[0,*,*]), reform(radec_img[1,*,*])) return, cxyz end ;mask keyword set if already read a fits mask, in order to save time ;function mask_cat_image, mask_file, image, radec_cent, ra_min, dec_min, $ ; spix, xstart, ystart, mask, proj_option = $ ; proj_option, astrometry = astrometry, test = test pro mask_cat_image, mask_file, image, radec_cent, ra_min, dec_min, $ spix, xstart, ystart, mask, mask_hdr, non_masked, $ proj_option = proj_option, astrometry = astrometry, $ test = test ;by default, catalog coordinates are assumed to be already projected ;at equator (if not, set keyword proj_option to 1) if(not keyword_set(proj_option)) then proj_option = 0 image_s = size(image, /dimensions) nxpix = image_s[0] nypix = image_s[1] ;check whether mask file is a .fits image check_fits_mask = strmatch(mask_file, '*fits*') ;mask is not a fits file, so assumes it is an ascii one if(check_fits_mask eq 0) then begin ;read mask file and mask pixels convex_true = 0 line = '' ilines = 0L openr, 44, mask_file while(not eof(44)) do begin readf, 44,line if (strmatch(line, 'CONVEX*') eq 1) then begin if(convex_true eq 0) then $ cxyz_img = get_cxyz_image(nxpix, nypix, ra_min, dec_min, spix, $ astrometry = astrometry, radec_cent = $ radec_cent) masked_line = mask_convex(line, cxyz_img) convex_true = 1 endif else begin if (strmatch(line, 'CIRCLE*') eq 0) then begin if (keyword_set(astrometry)) then begin masked_line = $ mask_poly(line, nxpix, nypix, ra_min, dec_min, spix, $ astrometry = astrometry) endif else begin masked_line = $ mask_poly(line, nxpix, nypix, ra_min, dec_min, spix, $ radec_cent = radec_cent) endelse endif else begin if (keyword_set(astrometry)) then begin masked_line = mask_circle(line, nxpix, nypix, ra_min, $ dec_min, spix, astrometry = $ astrometry, radec_cent = $ radec_cent, image = image) endif else begin if(keyword_set(test)) then begin masked_line = mask_circle(line, nxpix, nypix, ra_min, $ dec_min, spix, /test) endif else begin masked_line = mask_circle(line, nxpix, nypix, ra_min, $ dec_min, spix) endelse endelse endelse endelse ;deal with inclusion mask if (strmatch(line, '*_IN*') eq 1 ) then begin masked_line = setdifference(lindgen(nxpix*nypix), masked_line) print, 'mask_in' endif if (ilines eq 0L) then begin masked = masked_line endif else begin masked = [masked, masked_line] endelse ilines++ endwhile close, 44 inside_mask = intarr(nxpix*nypix) inside_mask[masked] = 1 non_masked = where(inside_mask eq 0) endif else begin ;fits mask if(n_elements(mask) eq 0) then begin mask = mrdfits(mask_file, 0, mask_hdr) nxpix_m = sxpar(mask_hdr, 'NAXIS1') nypix_m = sxpar(mask_hdr, 'NAXIS2') endif else begin nxpix_m = n_elements(mask[*,0]) nypix_m = n_elements(mask[0,*]) endelse ;print, 'in mask_cat_image' ;help, mask ;help, mask_hdr ntotpix_m = nxpix_m*nypix_m ;get mask WCS astrometry extast, mask_hdr, mask_astr ;select rectangle in mask image according to xystart parameters ;mask_use = select_rect_image_no_astr(mask, xstart, ystart) ;get x, y pixel array for input image x_arr = lindgen(nxpix) y_arr = lindgen(nypix) xpix = cmreplicate(x_arr, [nypix]) ypix = transpose(cmreplicate(reverse(y_arr), [nxpix])) ;get ra, dec coordinates from input cat xy2ad, xpix, ypix, astrometry, ra, dec ;get x, y pixels in mask frame ad2xy, ra, dec, mask_astr, xpix_m, ypix_m ;get one dimension subscripts ;xypix_m = xpix_m + ypix_m*nxpix_m ;keep only pixels within mask subscript range ;use = where(xypix_m GE 0 and xypix_m LT ntotpix_m) use_xm = where(xpix_m GE 0 and xpix_m LT nxpix_m) use_ym = where(ypix_m[use_xm] GE 0 and ypix_m[use_xm] LT nypix_m) ;non_masked_temp = where(mask[xpix_m[use_xm], ypix_m[use_xm[use_ym]]] eq 0) non_masked_temp = where(mask[xpix_m[use_xm[use_ym]], $ ypix_m[use_xm[use_ym]]] eq 0) ;get one dimension subscripts for input image pixels non_masked_tx = use_xm[non_masked_temp] mod nxpix non_masked_ty = nypix-use_xm[use_ym[non_masked_temp]]/nxpix ;here nypix-... because fits y=0 pixel ;is on the bottom of the image, while ;IDL's is on the top non_masked = non_masked_tx+nxpix*non_masked_ty ;stop endelse ;return, non_masked end