pro get_astr_st, spix, ra_cent, dec_cent, x_cent, y_cent, astr

tagnames    = ['cd', 'cdelt', 'crpix', 'crval', 'ctype', 'longpole', $
               'latpole','pv2']
tagdescript = 'd(2,2),d(2),f(2),d(2),a(2),d,d,d(2)'
create_struct, astr, '', tagnames, tagdescript

astr.cd[0,0]  = 1.
astr.cd[1,1]  = 1.
astr.cd[1,0]  = -0.
astr.cdelt[0] = -spix
astr.cdelt[1] = spix
astr.crpix[0] = x_cent
astr.crpix[1] = y_cent
astr.crval[0] = ra_cent
astr.crval[1] = dec_cent
astr.ctype[0] = 'RA---TAN'
astr.ctype[1] = 'DEC--TAN'
astr.longpole = 180.
astr.latpole  = 90.
astr.pv2      = [0., 0.]

end

pro fill_grid_cat_xwtheta, infile, radec_cent, mask_file, rmin,  rmax, $
                           nbin, img, npix_ra, npix_dec, ra_min, dec_min, $
                           ra_max, dec_max, spix, xstart, ystart, custom_res,$
                           non_masked_img

;read cat file
cat = read_file(infile)

;project coordinates at equator
radec = radec2deg(radec_cent[0], radec_cent[1], cat[0,*], cat[1,*])


img = dblarr(npix_ra, npix_dec)

;compute pix coordinates for each object
xpix = long((radec[0,*] - ra_min)/spix)
ypix = long((radec[1,*] - dec_min)/spix)

;define a pseudo astrometry structure to handle masking
;coordinates of center in pixels
x_cent = long(-ra_min/spix)
y_cent = long(-dec_min/spix)
get_astr_st, spix, radec_cent[0], radec_cent[1], x_cent, y_cent, astr



cell_number = xpix+ ypix*npix_ra

ncells = npix_ra*npix_dec


;get the number of objects in each cell
histo_cells  = histogram(cell_number, max = ncells-1L, min = 0L)

;fill the image
img[lindgen(ncells)] = histo_cells


non_masked_img = mask_cat_image(mask_file, img, radec_cent, ra_min, $
                                dec_min, spix, xstart, ystart, astrometry = $
                                astr)
img_temp = dblarr(npix_ra, npix_dec)
img_temp[non_masked_img] = img[non_masked_img]
img = img_temp
;stop
end