Image and cube resampling

A common problem in astronomy is the resampling of 2D images and 3D cubes onto a common grid. Ideally, this is done only once in the data reduction workflow as each sub-pixel resampling redistributes the flux and creates correlation.

Assuming that one would like to re-grid and stack dithered images/cubes with a pre-determined World Coordinate System (WCS), two scenarios are supported:

  • Combine all images/cubes in a single cloud of points (a table) and when creating the resampled output image/cube all information is used in a single interpolation step, i.e. no image-by-image interpolation is done.

  • Resample each individual image/cube on the same output grid and combine the single resampled images in a second step with HDRL stacking methods.

Both scenarios interpolate only once. The second scenario is important in case not all bad pixels could be properly determined and inserted in the bad pixel mask upfront.

In PyHDRL this is hdrl.func.Resample, together with hdrl.func.ResampleMethod and hdrl.func.ResampleOutgrid.

Interpolation methods

The implemented interpolation algorithms are based on the MUSE pipeline and work for 2D images and 3D cubes:

  • hdrl.func.ResampleMethod.Nearest(): nearest neighbour resampling

  • hdrl.func.ResampleMethod.Linear(loop_distance, use_errorweights): inverse distance weighting

  • hdrl.func.ResampleMethod.Quadratic(loop_distance, use_errorweights): quadratic inverse distance weighting

  • hdrl.func.ResampleMethod.Renka(loop_distance, use_errorweights, critical_radius): Renka weighting function

  • hdrl.func.ResampleMethod.Drizzle(loop_distance, use_errorweights, pix_frac_x, pix_frac_y, pix_frac_lambda): drizzle-like weighting

  • hdrl.func.ResampleMethod.Lanczos(loop_distance, use_errorweights, kernel_size): Lanczos-like restricted sinc

Output grid

hdrl.func.ResampleOutgrid defines the resampled image or cube:

  • User2D(delta_ra, delta_dec, ra_min, ra_max, dec_min, dec_max, fieldmargin) / User3D(delta_ra, delta_dec, delta_lambda, ra_min, ra_max, dec_min, dec_max, lambda_min, lambda_max, fieldmargin): user-specified bounds and step sizes

  • Auto2D(delta_ra, delta_dec) / Auto3D(delta_ra, delta_dec, delta_lambda): only the step sizes are required; the rest is derived from the data

Compute

The routine does not work directly on an image or cube but on a cpl.core.Table (restable). For 2D images the table is created by hdrl.func.Resample.image_to_table() from an hdrl.core.Image and a cpl.drs.WCS; for a 3D data cube hdrl.func.Resample.imagelist_to_table() is used with an hdrl.core.ImageList and a cpl.drs.WCS.

If many images or cubes have to be combined into a single mosaic, the two functions can be called multiple times and the returned tables should be merged into a single table using cpl.core.Table.insert().

If the pixel to sky mapping cannot be encoded by a cpl.drs.WCS object, hdrl.func.Resample.restable_template(nrows) creates a table template to be filled by the pipeline developer.

compute() returns a hdrl.func.ResampleResult with hdr (cpl.core.PropertyList) and imlist (hdrl.core.ImageList).

wcs = cpl.drs.WCS(plist)
table = hdrl.func.Resample.image_to_table(himg, wcs)
rmethod = hdrl.func.ResampleMethod.Lanczos(1, False, 2)
outgrid = hdrl.func.ResampleOutgrid.Auto2D(0.01, 0.01)
result = hdrl.func.Resample.compute(table, rmethod, outgrid, wcs)
first_img = result.imlist[0].image
first_err = result.imlist[0].error
hdr = result.hdr