;If keywords ra_cent or dec_cent not set, look in the header of the
;image for the coordinates of the center of the field
;Needs hdr with WCS infos ...
;native_res: if set, use image native resolution, and do not resample
;image according to the binning
pro fill_grid_image, fits_file, mask_file, rmin, rmax, nbin, img, $
                     ra_min, dec_min, ra_max, dec_max, spix, non_masked_img, $
                     proj_corner, ra_cent = ra_cent, dec_cent = dec_cent, $
                     threshold = threshold, native_res = native_res, $
                     xstart = xstart, ystart = ystart, custom_res = custom_res
im_in = mrdfits(fits_file, 0, hdr)

if(not keyword_set(xstart)) then xstart = 0L
if(not keyword_set(ystart)) then ystart = 0L


;Get wcs infos
extast,hdr,astr
;Get ra and dec of the center if not supplied
if (not keyword_set(ra_cent)) then begin
    ra_cent  = sxpar(hdr, 'avaspra')
    dec_cent = sxpar(hdr, 'avaspdec')
endif

;select rectangle in image according to xystart parameters, and update
;WCS parameters
select_rect_image, im_in, xstart, ystart, im_use, astr

;First, mask image according to mask_file
non_masked_im = mask_cat_image(mask_file, im_use, [ra_cent, dec_cent], $
                               ra_min, dec_min,spix, xstart, ystart, $
                               astrometry = astr)



im_size = size(im_use, /dimensions)
nxpix = im_size[0]
nypix = im_size[1]

;create masked image
img_masked                = dblarr(nxpix, nypix)
img_masked[non_masked_im] = im_use[non_masked_im]



;create mask image
img_mask                = intarr(nxpix, nypix)
img_mask[non_masked_im] = 1

non_masked_im_cols = non_masked_im mod nxpix
non_masked_im_rows = non_masked_im/nxpix


;get (ra,dec) coords of non masked pixels
xy2ad, non_masked_im_cols, non_masked_im_rows, astr, non_masked_ra, $
  non_masked_dec

;project coordinates at equator
coord_proj = radec2deg(ra_cent, dec_cent, non_masked_ra, non_masked_dec)

 ;Define grid
ra_max   = max(coord_proj[0,*], min = ra_min)
dec_max  = max(coord_proj[1,*], min = dec_min)
ra_max  += 0.075D
dec_max += 0.075D
ra_min  -= 0.075D
dec_min -= 0.075D



delta_ra = ra_max-ra_min
delta_dec = dec_max-dec_min
    
rmin_log  = alog10(rmin)
rmax_log  = alog10(rmax)
bin_log   = (rmax_log-rmin_log)/nbin
r_arr     = 10^(rmin_log +bin_log*dindgen(nbin))
dr        = r_arr[1] - r_arr[0]


  
;handle resolution

if(keyword_set(native_res)) then begin
;native resolution
   ;use input image size
   npix_ra  = nxpix
   npix_dec = nypix
   spix     = astr.cdelt[1]
endif else if(keyword_set(custom_res)) then begin
   ;use custom resolution factor
   npix_ra  = nxpix*1./custom_res
   npix_dec = nypix*1./custom_res
   spix     = astr.cdelt[1]*custom_res
endif else begin
    ;use binning to define image resampling
   npix_ra  = long(delta_ra/dr) + 1L
   npix_dec = long(delta_dec/dr) + 1L
   spix     = delta_ra/npix_ra
endelse


;resample image and image mask if needed
if(not keyword_set(native_res)) then begin
   img  = congrid(img_masked, npix_ra, npix_dec)
   mask = congrid(img_mask, npix_ra, npix_dec)
endif else begin
   img  = img_masked
   mask = img_mask
endelse

non_masked_img = where(mask ne 0)

;apply threshold if needed
if (keyword_set(threshold)) then begin
   non_masked_img = apply_img_threshold(img, non_masked_img, threshold)
endif


;get coordinates of bottom left corner
xy2ad, 0, 0, astr, corner_ra, corner_dec
;projected coordinates of bottom left corner
proj_corner = radec2deg(ra_cent, dec_cent, corner_ra, corner_dec)
end