Library demo - Warp to another image (GNU Astronomy Utilities)

Next: Library demo - Warp to new grid, Previous: Library demo - reading and writing table columns, Up: Library demo programs [Contents][Index]

12.4.5 Library demo - Warp to another image ¶

Gnuastro’s warp library (that you can access by including gnuastro/warp.h) allows you to resample an image from a grid to another entirely using the WCSLIB (while accounting for distortions if necessary; see Warp library (warp.h)). The Warp library uses a pixel-mixing or area-based resampling approach which is fully described in Resampling. The most generic uses cases for this library are already available in the Invoking Warp program. For a related demo (where the output grid and WCS are constructed from scratch), see Library demo - Warp to new grid.

In the example below, we are warping the input.fits file to the same pixel grid and WCS as reference.fits image (assuming it is in hdu 0). You can download the FITS files in the Color channels in same pixel grid section and use them as input.fits and reference.fits files. Feel free to change these names to your own test file names. This can be useful when you have a complex grid and WCS containing various keywords such as non-linear distortion coefficients, etc. For example datasets, see the description of the --gridfile option in Align pixels with WCS considering distortions.

To compile the demonstration program below, copy and paste the contents in a plain-text file (let’s assume you named it align-to-img.c) and use BuildProgram with this command: ‘astbuildprog align-to-img.c’. Please note that the demo program does not perform many sanity checks to avoid making it too complex and to highlight this particular feature in the library. For a robust method write programs with all the necessary sanity checks, see Gnuastro’s Warp source code, see Program source.

To encourage good coding practices, this script contains a copyright notice with a place holder for your name and your email (as you customize it for your own purpose). Always keep a one-line description and copyright notice like this in all your scripts, such “metadata” is very important to accompany every source file you write. Of course, when you write the source file from scratch and just learn how to use a single function from this manual, only your name/year should appear. The existing name of the original author of this example program is only for cases where you copy-paste this whole file.

/* Warp to another image.
 *
 * In the example below, we are warping the input.fits file to the same
 * pixel grid and WCS as reference.fits image.
 *
 * Copyright (C) 2024      Your Name <your@@email.address>
 * Copyright (C) 2022-2024 Pedram Ashofteh-Ardakani <pedramardakani@pm.me>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>

#include <gnuastro/wcs.h>       /* contains gnuastro's fits.h */
#include <gnuastro/warp.h>      /* contains gnuastro's data.h */
#include <gnuastro/array.h>     /* contains gnuastro's type.h */

int
main(void)
{
  /* Input file's name and HDU. */
  char *filename="input.fits", *hdu="1";

  /* Reference file's name and HDU. */
  char *gridfile="reference.fits", *gridhdu="0";

  /* Output file name. */
  char *outname="align-to-img.fits";

  /* Low-level variables needed to read the reference file's size. */
  int nwcs;
  size_t ndim, *dsize;

  /* Initialize the 'wa' struct with empty values and NULL pointers. */
  gal_warp_wcsalign_t wa=gal_warp_wcsalign_template();

  /* Read the input image and its WCS. */
  wa.input=gal_array_read_one_ch_to_type(filename, hdu, NULL,
                                         GAL_TYPE_FLOAT64, -1,  0, NULL);
  wa.input->wcs=gal_wcs_read(filename, hdu, 0, 0, 0, &wa.input->nwcs,
                             NULL);

  /* Prepare the warp input structure, use all threads available. */
  wa.coveredfrac=1; wa.edgesampling=0; wa.numthreads=0;

  /* Set the target grid to be the same as wcsref.fits file on hdu 0. */
  wa.twcs=gal_wcs_read(gridfile, gridhdu, 0, 0, 0, &nwcs, NULL);
  if(wa.twcs==NULL)
    {
      fprintf(stderr, "%s (hdu %s): no WCS! Can't continue\n",
	      gridfile, gridhdu);
      exit(EXIT_FAILURE);
    }

  /* Read the output image size (from the reference image). Note that
   * 'dsize' will be freed while freeing 'widthinpix'). */
  dsize=gal_fits_img_info_dim(gridfile, gridhdu, &ndim, NULL);

  /* Convert the 'dsize' to a 'gal_data_t' so the library can use it. */
  wa.widthinpix=gal_data_alloc(dsize, GAL_TYPE_SIZE_T, 1, &ndim,
			       NULL, 1, -1, 0, NULL, NULL, NULL);

  /* Do the warp, then convert the output to a 32-bit float (the default
   * float64 is too much for observational data and just wastes
   * storage!). But if you are warping mock data before adding noise
   * (where you do have float64 level precision), remove the type
   * conversion line. */
  gal_warp_wcsalign(&wa);
  wa.output=gal_data_copy_to_new_type_free(wa.output, GAL_TYPE_FLOAT32);

  /* WARNING: make sure there is no file with same name as 'out.fits'
   * or the result will be appended to its final HDU. */
  gal_fits_img_write(wa.output, outname, NULL, 0);

  /* Clean up. */
  gal_data_free(wa.input);
  gal_data_free(wa.output);
  gal_data_free(wa.widthinpix);

  /* Give control back to the operating system. */
  return EXIT_SUCCESS;
}