[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 70693618: Warp: aligning the pixel and WCS coo
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 70693618: Warp: aligning the pixel and WCS coords, while correcting distortions |
|
Date: |
Sun, 2 Oct 2022 20:24:47 -0400 (EDT) |
branch: master
commit 70693618e9dfc8eb1fdb8ce04e53a2179326fdcf
Author: Pedram Ashofteh Ardakani <pedramardakani@pm.me>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Warp: aligning the pixel and WCS coords, while correcting distortions
Until now, Warp could only do linear resampling of the input grid. However,
astronomical imaging cameras often have distortions that also need to be
corrected prior to stacking of the images into one pixel grid. Furthermore
the '--align' operator (which tried to estimate the necessary linear
rotation on images with no distortion) would also fail in some peculiar
scenarios.
With this commit, Warp will now construct the output pixel grid and WCS
based on the user desired properties and correct any form of distortion in
the process. This is now the default operation of Warp (when no linear
operations are requested). Until now, Warp would crash with an error if no
linear operations were given.
Unlike Warp's previous linear warps, the new distortion correcting aligning
operations of warp have been written in a newly added library in Gnuastro:
'warp.h'. The Warp program calls these library functions for this job. But
as a library, it can be used in many other scenarios also, not just the
Warp program.
Furthermore, given that the Warp program was one of the first written
programs in Gnuastro, some aspects of it weren't "modern" (for example it
didn't use the high-level 'gnuastro/threads.h' library features and had a
low-level implementation of it, or it didn't follow the more recent naming
conventions in its function names). Such things have also been corrected.
The following changes have also been made with this commit:
- The linear '--align' option has been removed: this is now the default
behavior of Warp.
- Some memory leaks have been found and corrected.
---
NEWS | 39 +-
bin/warp/args.h | 89 ++++-
bin/warp/astwarp.conf | 6 +-
bin/warp/main.c | 2 +-
bin/warp/main.h | 11 +-
bin/warp/ui.c | 489 +++++++++++++-----------
bin/warp/ui.h | 17 +-
bin/warp/warp.c | 286 +++++++-------
doc/gnuastro.texi | 878 ++++++++++++++++++++++++++++++++----------
lib/Makefile.am | 3 +
lib/gnuastro/dimension.h | 23 ++
lib/gnuastro/warp.h | 109 ++++++
lib/gnuastro/wcs.h | 3 +
lib/warp.c | 970 +++++++++++++++++++++++++++++++++++++++++++++++
lib/wcs.c | 48 ++-
tests/during-dev.sh | 2 -
16 files changed, 2354 insertions(+), 621 deletions(-)
diff --git a/NEWS b/NEWS
index bb466aba..0f8609bc 100644
--- a/NEWS
+++ b/NEWS
@@ -104,6 +104,20 @@ See the end of the file for license conditions.
--fitestimatehdu: HDU containing table in file given to '--fitestimate'.
--fitestimatecol: Column containing X axis values for '--fitestimate'.
+ Warp:
+ - Can correct distortions (with any standard recognized by WCSLIB) and
+ simultaneously align the image to the coordinate system. When no named
+ linear operation (like '--rotate', '--scale' or etc) is requested, Warp
+ will go into this mode. It is highly customizable through the following
+ options. See the "Invoking Warp" section of the book for more.
+ --center: RA, DEC of the center of the central pixel of output.
+ --widthinpix: Width of output in pixels.
+ --cdelt: Pixel scale of output ('CDELTi' keywords in FITS).
+ --ctype: Coordinates and projection algorithm. Default: RA/Dec and
+ Gnomonic or 'TAN').
+ --edgesampling: extra sampling of pixel polygon to account for strong
+ non-linear projection or distortion effects, when necessary.
+
astscript-fits-view:
--ds9colorbarmulti: show a separate color-bar for each image in DS9. By
default this script will show a single color-bar for all the images to
@@ -188,15 +202,26 @@ See the end of the file for license conditions.
- gal_units_mag_to_sb: surface brightness (SB) from magnitude and area.
- gal_units_sb_to_counts: counts from SB, zeropoint and area.
- gal_units_sb_to_mag: magnitude from SB and area.
+ - gal_warp_wcsalign_init: initialize the WCS aligning structure.
+ - gal_warp_wcsalign_onpix: Per-pixel filling of output.
+ - gal_warp_wcsalign_onthread: function to give to pthreads.
+ - gal_warp_wcsalign: high-level function to align input by its WCS.
+ - gal_warp_wcsalign_free: free the contents of the WCS aligning structure.
+ - gal_wcs_free: free a WCS structure that is created or read by Gnuastro.
** Removed features
Statistics:
- --refcol has been removed because it breaks the modularity principle
- (given that it is the job of Gnuastro's Table program to limit rows
- from a larger table based on many different criteria). The output of
- Table can be directly piped to Statistics to achieve the same (and
- much more feature-rich effect).
+ --refcol has been removed because it breaks the modularity principle
+ (given that it is the job of Gnuastro's Table program to limit rows
+ from a larger table based on many different criteria). The output of
+ Table can be directly piped to Statistics to achieve the same (and much
+ more feature-rich effect).
+
+ Warp:
+ --align: has been removed. This is because aligning an image (while
+ correcting for any possible distortion) is now the default behavior of
+ Warp (when no linear operations have been requested).
** Changed features
@@ -206,6 +231,10 @@ See the end of the file for license conditions.
level. The "Sufi simulates a detection" (which was previously first)
has been moved to the fourth section.
+ Warp:
+ - The short format of the '--centeroncorner' option has been removed. The
+ '-c' is now the short format for the new '--center' option to Warp.
+
astscript-psf-scale-factor:
--widthinpix: new name for the old '--stampwidth' option. This was done
to have the same name to a similar option in Crop and help in
diff --git a/bin/warp/args.h b/bin/warp/args.h
index 841b81a0..5d476ea4 100644
--- a/bin/warp/args.h
+++ b/bin/warp/args.h
@@ -5,6 +5,7 @@ Warp is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <mohammad@akhlaghi.org>
Contributing author(s):
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
Copyright (C) 2016-2022 Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
@@ -62,6 +63,8 @@ struct argp_option program_options[] =
+
+
/* Output. */
{
"keepwcs",
@@ -86,29 +89,97 @@ struct argp_option program_options[] =
&p->coveredfrac,
GAL_TYPE_FLOAT64,
GAL_OPTIONS_RANGE_GE_0_LE_1,
- GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+
+
+
{
0, 0, 0, 0,
- "Warps:",
- UI_GROUP_WARPS
+ "Align with WCS coordinates (correcting distortion, default mode)",
+ UI_GROUP_ALIGN
},
{
- "align",
- UI_KEY_ALIGN,
+ "center",
+ UI_KEY_CENTER,
+ "FLT,FLT",
0,
+ "Center coordinate of the output image in RA,DEC.",
+ UI_GROUP_ALIGN,
+ &p->wa.center,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_parse_csv_float64
+ },
+ {
+ "widthinpix",
+ UI_KEY_WIDTHINPIX,
+ "INT,INT",
0,
- "Align the image and celestial axes.",
- UI_GROUP_WARPS,
+ "Output image width and height in pixels.",
+ UI_GROUP_ALIGN,
+ &p->wa.widthinpix,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_GT_0_ODD,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_parse_csv_float64
+ },
+ {
+ "ctype",
+ UI_KEY_CTYPE,
+ "STR[,STR]",
0,
- GAL_TYPE_INVALID,
+ "FITS standard CTYPE value (e.g., 'RA---TAN').",
+ UI_GROUP_ALIGN,
+ &p->wa.ctype,
+ GAL_TYPE_STRING,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET,
- ui_add_to_modular_warps_ll
+ gal_options_parse_csv_strings
+ },
+ {
+ "cdelt",
+ UI_KEY_CDELT,
+ "FLT[,FLT]",
+ 0,
+ "Pixel scale of output (usually degrees/pixel).",
+ UI_GROUP_ALIGN,
+ &p->cdelt,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_GT_0,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_parse_csv_float64
+ },
+ {
+ "edgesampling",
+ UI_KEY_EDGESAMPLING,
+ "INT",
+ 0,
+ "Number of extra samplings in pixel sides.",
+ UI_GROUP_ALIGN,
+ &p->wa.edgesampling,
+ GAL_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GE_0,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ },
+
+
+
+
+
+ {
+ 0, 0, 0, 0,
+ "Linear warps (must be called explicitly on command-line)",
+ UI_GROUP_WARPS
},
{
"rotate",
diff --git a/bin/warp/astwarp.conf b/bin/warp/astwarp.conf
index 549e7286..86c004b0 100644
--- a/bin/warp/astwarp.conf
+++ b/bin/warp/astwarp.conf
@@ -20,9 +20,11 @@
# warranty.
# Input:
+ edgesampling 0
# Output:
- coveredfrac 1.0
+ coveredfrac 1.0
+ ctype RA---TAN,DEC--TAN
# Common parameters
- type float32
+ type float32
diff --git a/bin/warp/main.c b/bin/warp/main.c
index 0d52a922..c63047fe 100644
--- a/bin/warp/main.c
+++ b/bin/warp/main.c
@@ -36,7 +36,7 @@ int
main (int argc, char *argv[])
{
struct timeval t1;
- struct warpparams p={{{0},0},0};
+ struct warpparams p={{{0},0},{0},0};
/* Set the starting time.*/
time(&p.rawtime);
diff --git a/bin/warp/main.h b/bin/warp/main.h
index fc6f25df..c0ba6e0b 100644
--- a/bin/warp/main.h
+++ b/bin/warp/main.h
@@ -5,6 +5,7 @@ Warp is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <mohammad@akhlaghi.org>
Contributing author(s):
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
Copyright (C) 2016-2022 Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
@@ -25,6 +26,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Include necessary headers */
#include <gnuastro/data.h>
+#include <gnuastro/warp.h>
#include <gnuastro-internal/options.h>
@@ -42,26 +44,29 @@ struct warpparams
{
/* From command-line */
struct gal_options_common_params cp; /* Common parameters. */
+ gal_warp_wcsalign_t wa; /* Nonlinear-specific parameters. */
char *inputname; /* Name of input file. */
size_t hstartwcs; /* Header keyword No. to start reading WCS. */
size_t hendwcs; /* Header keyword No. to end reading WCS. */
uint8_t keepwcs; /* Wrap the warped/transfomed pixels. */
uint8_t centeroncorner; /* Shift center by 0.5 before and after. */
double coveredfrac; /* Acceptable fraction of output covered. */
+ gal_data_t *cdelt; /* Pixel scale of the output image. */
/* Internal parameters: */
gal_data_t *input; /* Input data structure. */
gal_data_t *output; /* output data structure. */
gal_data_t *matrix; /* Warp/Transformation matrix. */
gal_data_t *modularll; /* List of modular warpings. */
- double *inverse; /* Inverse of the input matrix. */
double *inwcsmatrix; /* Input WCS matrix. */
- double *pixelscale; /* Pixel scale of input image. */
+ double *inverse; /* Inverse of the input matrix. */
time_t rawtime; /* Starting time of the program. */
- size_t extinds[4]; /* Indexs of the minimum and maximum values. */
size_t ordinds[4]; /* Indexs of anticlockwise vertices. */
+ size_t extinds[4]; /* Indexs of the minimum and maximum values. */
double outfpixval[2]; /* Pixel value of first output pixel. */
double opixarea; /* Area of output pix in units of input pix. */
+ uint8_t nonlinearmode; /* If warp must work in nonlinear mode. */
+ uint8_t distortiontype; /* Store distortion type in nonlinear mode. */
};
#endif
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index de62e4e6..e88ed0f8 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -5,6 +5,7 @@ Warp is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <mohammad@akhlaghi.org>
Contributing author(s):
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
Copyright (C) 2016-2022 Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
@@ -29,6 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <gnuastro/wcs.h>
+#include <gnuastro/warp.h>
#include <gnuastro/fits.h>
#include <gnuastro/array.h>
#include <gnuastro/table.h>
@@ -65,10 +67,12 @@ static char
args_doc[] = "ASTRdata";
const char
-doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will warp/transform the "
- "input image using an input coordinate matrix. Currently it accepts any "
- "general projective mapping (which includes affine mappings as a "
- "subset). \n"
+doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will resample the pixel "
+ "grid of an input image. By default (if no special linear warping is "
+ "requested), it will align the image to the WCS coordinates in any"
+ "and remove any possible distortion. Linear warping (like '--rotate' "
+ "or '--scale') should be explicitly requested with the options under "
+ "the \"Linear warps\" group below. \n"
GAL_STRINGS_MORE_HELP_INFO
/* After the list of options: */
"\v"
@@ -115,6 +119,8 @@ ui_initialize_options(struct warpparams *p,
cp->numthreads = gal_threads_number();
cp->coptions = gal_commonopts_options;
+ /* Program specific initializations. */
+ p->wa.edgesampling = GAL_BLANK_SIZE_T;
/* Set the mandatory common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -127,6 +133,7 @@ ui_initialize_options(struct warpparams *p,
break;
case GAL_OPTIONS_KEY_SEARCHIN:
+ case GAL_OPTIONS_KEY_IGNORECASE:
case GAL_OPTIONS_KEY_TABLEFORMAT:
case GAL_OPTIONS_KEY_STDINTIMEOUT:
cp->coptions[i].flags=OPTION_HIDDEN;
@@ -225,34 +232,15 @@ ui_add_to_modular_warps_ll(struct argp_option *option,
char *arg,
gal_data_t *new;
struct warpparams *p=(struct warpparams *)params;
- /* When an argument is necessary (note that '--align' doesn't take an
- argument), make sure we actually have a string to parse. Also note
- that if an argument is necessary, but none is given Argp will
- automatically abort the program with an informative error. */
+ /* When an argument is necessary, make sure we actually have a string to
+ parse. Also note that if an argument is necessary, but none is given
+ Argp will automatically abort the program with an informative
+ error. */
if(arg && *arg=='\0')
error(EXIT_FAILURE, 0, "empty string given to '--%s'", option->name);
/* Parse the (possible) arguments. */
- if(option->key==UI_KEY_ALIGN)
- {
- /* For functions the standard checking isn't done, so first, we'll
- make sure that if we are in a configuration file (where
- 'arg!=NULL'), the value is either 0 or 1. */
- if( arg && strcmp(arg, "0") && strcmp(arg, "1") )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "the '--align' "
- "option takes no arguments. In a configuration file "
- "it can only have the values '1' or '0', indicating "
- "if it should be used or not");
-
- /* Align doesn't take any values, but if called in a configuration
- file with a value of '0', we should ignore it. */
- if(arg && *arg=='0') return NULL;
-
- /* Allocate the data structure. */
- new=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 0, NULL, NULL, 0, -1, 1,
- NULL, NULL, NULL);
- }
- else new=gal_options_parse_list_of_numbers(arg, filename, lineno);
+ new=gal_options_parse_list_of_numbers(arg, filename, lineno);
/* If this was a matrix, then put it in the matrix element of the main
@@ -334,53 +322,231 @@ ui_add_to_modular_warps_ll(struct argp_option *option,
char *arg,
/*************** Sanity Check *******************/
/**************************************************************/
static void
-ui_check_options_and_arguments(struct warpparams *p)
+ui_check_options_and_arguments_nonlinear(struct warpparams *p)
{
- /* Read the input.*/
- if(p->inputname)
+ gal_warp_wcsalign_t *wa=&p->wa;
+ size_t two=2, *sarray=NULL, stemp, indim=0;
+ double *icenter=NULL, *iwidth=NULL, *imin=NULL, *imax=NULL, *tmp=NULL;
+
+ /* Put the number of threads into the non-linear structure. */
+ wa->numthreads=p->cp.numthreads;
+
+ /* Sanity check user's possibly given '--center'. */
+ if(wa->center)
{
- /* Make sure a HDU is given. */
- if( gal_fits_file_recognized(p->inputname) && p->cp.hdu==NULL )
- error(EXIT_FAILURE, 0, "no HDU specified, you can use the '--hdu' "
- "('-h') option and give it the HDU number (starting from "
- "zero), or extension name (generally, anything acceptable "
- "by CFITSIO)");
-
- /* Read the input image as double type and its WCS structure. */
- p->input=gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu,
- NULL, GAL_TYPE_FLOAT64,
- p->cp.minmapsize,
- p->cp.quietmmap);
-
- /* Read the WCS and remove one-element wide dimension(s). */
- p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu,
- p->cp.wcslinearmatrix, p->hstartwcs,
- p->hendwcs, &p->input->nwcs);
- p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
- p->input->dsize,
- p->input->wcs);
-
- /* Currently Warp only works on 2D images. */
- if(p->input->ndim!=2)
- error(EXIT_FAILURE, 0, "input has %zu dimensions but Warp currently "
- "only works on 2D datasets (images).\n\n"
- "We do plan to add 3D functionality (see "
- "https://savannah.gnu.org/task/?15729), so please get in "
- "touch if you need it (any further interest, support or help "
- "would be useful)", p->input->ndim);
-
- /* Get basic WCS information. */
- if(p->input->wcs)
+ if(wa->center->size !=2)
+ error(EXIT_FAILURE, 0, "%zu value(s) given to '--center', "
+ "however this option takes exactly 2 values to specify "
+ "the output image center", wa->center->size);
+
+ tmp=wa->center->array;
+ if(tmp[0] < 0 || tmp[0] > 360 || tmp[1] < -90 || tmp[1] > 90 )
+ error(EXIT_FAILURE, 0, "the first value '--center' should be "
+ "between 0 and 360 (inclusive, because it is the RA) "
+ "and the second value should be between -90 and 90 "
+ "(inclusive, because it is the Dec), however the given "
+ "values are: %g and %g", tmp[0], tmp[1]);
+ }
+ else /* --center not given, use the input image to define it. */
+ {
+ /* Get sky coverage information about the input image. Note that
+ this allocates the pointers that have to be freed later. */
+ if( gal_wcs_coverage(p->inputname, p->cp.hdu, &indim, &icenter,
+ &iwidth, &imin, &imax)==0 )
+ error(EXIT_FAILURE, 0, "%s (hdu %s): is not usable for finding "
+ "sky coverage", p->inputname, p->cp.hdu);
+
+
+ /* Store the center array and clean up the rest. */
+ wa->center=gal_data_alloc(icenter, GAL_TYPE_FLOAT64, 1,
+ &two, NULL, 0, p->cp.minmapsize,
+ p->cp.quietmmap, NULL, NULL,
+ NULL);
+
+ free(imin);
+ free(imax);
+ free(iwidth);
+ }
+
+ /* If the output width is given, make sure it has exactly two
+ values for width and height. */
+ if(wa->widthinpix)
+ {
+ if( wa->widthinpix->size != 2 )
+ error(EXIT_FAILURE, 0, "%zu value(s) given to '--widthinpix', "
+ "however this option takes exactly 2 values as the "
+ "output image width and height in pixels",
+ wa->widthinpix->size);
+
+ /* Image size must be of type size_t, but first we should make sure
+ the user didn't give a floating-point value. */
+ tmp=wa->widthinpix->array;
+ if( tmp[0]!=ceil(tmp[0]) || tmp[1]!=ceil(tmp[1]) )
+ error(EXIT_FAILURE, 0, "value to '--widthinpix' must be "
+ "integers, but they are: %g, %g", tmp[0], tmp[1]);
+ wa->widthinpix=gal_data_copy_to_new_type(wa->widthinpix,
+ GAL_TYPE_SIZE_T);
+
+ /* Image size must be ODD */
+ sarray=wa->widthinpix->array;
+ if( sarray[0]%2==0 || sarray[1]%2==0 )
{
- p->pixelscale=gal_wcs_pixel_scale(p->input->wcs);
- if(p->pixelscale==NULL)
- error(EXIT_FAILURE, 0, "%s (hdu %s): the pixel scale couldn't "
- "be deduced from the WCS.", p->inputname, p->cp.hdu);
- p->inwcsmatrix=gal_wcs_warp_matrix(p->input->wcs);
+ /* Let the user know that we are updating the output size. */
+ error(EXIT_SUCCESS, 0, "WARNING: '--width' must be odd: "
+ "updating %zux%zu to %zux%zu", sarray[0], sarray[1],
+ sarray[0]%2 ? sarray[0] : sarray[0]+1,
+ sarray[1]%2 ? sarray[1] : sarray[1]+1);
+
+ /* Keep the final values. */
+ if(sarray[0]%2==0) sarray[0]+=1;
+ if(sarray[1]%2==0) sarray[1]+=1;
}
+
+ /* To keep the later steps consistent with the C ordering of
+ dimensions, we'll swap the fast and slow axis from FITS (which has
+ the same convention as FORTRAN), because the user sees a FITS
+ file. In an aligned image, sarray[0] is the horizontal axis and
+ sarray[1] the vertical one.
+
+ NAXIS1 -> dsize[1] -> horizontal axis
+ NAXIS2 -> dsize[0] -> vertical axis */
+ stemp=sarray[0]; sarray[0]=sarray[1]; sarray[1]=stemp;
+ }
+
+ /* Check CTYPE */
+ if(!wa->ctype)
+ error(EXIT_FAILURE, 0, "no output projection CTYPE specified, "
+ "you can use the '--ctype' option and give it a comma "
+ "separated CTYPE value recognized by the WCSLIB (e.g. "
+ "--ctype=RA---TAN,DEC--TAN)");
+ if(wa->ctype->size != p->input->ndim)
+ error(EXIT_FAILURE, 0, "%zu value(s) given to '--ctype', but it "
+ "takes exactly 2 values", wa->ctype->size);
+
+
+ /* Copy necessary parameters for the nonlinear warp to work */
+ wa->input=p->input;
+ wa->cdelt=p->cdelt;
+ wa->coveredfrac=p->coveredfrac;
+}
+
+
+
+
+static void
+ui_check_cdelt(struct warpparams *p)
+{
+ size_t two=2, i;
+ gal_data_t *olddata=NULL;
+ double *tmp=NULL, *cdelt=NULL;
+
+ if(!p->cdelt)
+ {
+ /* CDELT is not given, try to deduce from WCS */
+ cdelt=gal_wcs_pixel_scale(p->input->wcs);
+ if(!cdelt)
+ error(EXIT_FAILURE, 0, "%s (hdu %s): the pixel scale couldn't "
+ "be deduced from the WCS.", p->inputname, p->cp.hdu);
+
+ /* Set CDELT to the maximum value */
+ cdelt[0] = ( cdelt[0] > cdelt[1] ? cdelt[0] : cdelt[1] );
+ cdelt[1] = cdelt[0];
+ p->cdelt=gal_data_alloc(cdelt, GAL_TYPE_FLOAT64, 1, &two, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap, NULL, NULL,
+ NULL);
}
else
+ {
+ /* CDELT is given, make sure there are no more than two values */
+ if(p->cdelt->size > 2)
+ error(EXIT_FAILURE, 0, "%zu values given to '--cdelt', "
+ "however this option takes no more than 2 values",
+ p->cdelt->size);
+
+ /* If only one value given to CDELT */
+ if(p->cdelt->size == 1)
+ {
+ olddata=p->cdelt; tmp=olddata->array;
+ p->cdelt=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1,
+ &two, NULL, 0,
+ p->cp.minmapsize,
+ p->cp.quietmmap, NULL, NULL,
+ NULL);
+ cdelt=p->cdelt->array;
+ cdelt[0]=cdelt[1]=tmp[0];
+ gal_data_free(olddata);
+ }
+
+ /* Check if the CDELT is in a reasonable range of degrees */
+ cdelt=p->cdelt->array;
+ for(i=p->cdelt->size; i--;)
+ if(cdelt[i]>0.01)
+ error(EXIT_SUCCESS, 0, "WARNING: CDELT on dimension %zu has "
+ "the unusual value of %f degrees. If you meant to "
+ "define CDELT in arcseconds please use: "
+ "'--cdelt=%g/3600'", i, cdelt[i], cdelt[i]);
+ }
+}
+
+
+
+
+
+static void
+ui_check_options_and_arguments(struct warpparams *p)
+{
+ /* Read the input.*/
+ if(p->inputname==NULL)
error(EXIT_FAILURE, 0, "no input file is specified");
+
+ /* Make sure a HDU is given. */
+ if( gal_fits_file_recognized(p->inputname) && p->cp.hdu==NULL )
+ error(EXIT_FAILURE, 0, "no HDU specified, you can use the '--hdu' "
+ "('-h') option and give it the HDU number (starting from "
+ "zero), or extension name (generally, anything acceptable "
+ "by CFITSIO)");
+
+ /* Make sure mandatory options are provided. */
+ if(p->modularll==NULL && p->matrix==NULL)
+ {
+ p->nonlinearmode=1;
+ if(p->wa.edgesampling==GAL_BLANK_SIZE_T)
+ error(EXIT_FAILURE, 0, "no '--edgesampling' provided");
+ }
+
+ /* Read the input image as double type and its WCS structure. */
+ p->input=gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu,
+ NULL, GAL_TYPE_FLOAT64,
+ p->cp.minmapsize,
+ p->cp.quietmmap);
+
+ /* Read the WCS and remove one-element wide dimension(s). */
+ p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu,
+ p->cp.wcslinearmatrix, p->hstartwcs,
+ p->hendwcs, &p->input->nwcs);
+ p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
+ p->input->dsize,
+ p->input->wcs);
+
+ /* Currently Warp only works on 2D images. */
+ if(p->input->ndim!=2)
+ error(EXIT_FAILURE, 0, "input has %zu dimensions but Warp currently "
+ "only works on 2D datasets (images).\n\n"
+ "We do plan to add 3D functionality (see "
+ "https://savannah.gnu.org/task/?15729), so please get in "
+ "touch if you need it (any further interest, support or help "
+ "would be useful)", p->input->ndim);
+
+ /* Get basic WCS information. */
+ if(p->input->wcs)
+ {
+ ui_check_cdelt(p);
+ p->inwcsmatrix=gal_wcs_warp_matrix(p->input->wcs);
+ }
+
+ /* Do all the distortion correction sanity-checks.*/
+ if(p->nonlinearmode)
+ ui_check_options_and_arguments_nonlinear(p);
}
@@ -469,133 +635,6 @@ ui_matrix_prepare_raw(struct warpparams *p)
-/* Set the matrix so the image is aligned with the axises. Note that
- WCSLIB automatically fills the CRPI */
-static void
-ui_matrix_make_align(struct warpparams *p, double *tmatrix)
-{
- double A, *w, *P, x[4];
-
- /* Make sure the input image had a WCS structure. */
- if(p->input->wcs==NULL)
- error(EXIT_FAILURE, 0, "%s (hdu: %s): no WCS information present, "
- "hence the '--align' option cannot be used", p->inputname,
- p->cp.hdu);
-
- /* Check if there is only two WCS axises: */
- if(p->input->wcs->naxis!=2)
- error(EXIT_FAILURE, 0, "the WCS structure of %s (hdu: %s) has %d "
- "axises. For the '--align' option to operate it must be 2",
- p->inputname, p->cp.hdu, p->input->wcs->naxis);
-
-
- /* For help in reading, use aliases for the WCS matrix and pixel scale.*/
- P=p->pixelscale;
- w=p->inwcsmatrix;
-
-
- /* Lets call the given WCS orientation 'W', the rotation matrix we want
- to find as 'X' and the final (aligned matrix) 'P' (which is the pixel
- scale):
-
- X W P
- ------ ------ -------
- x0 x1 w0 w1 -P0 0
- x2 x3 * w2 w3 = 0 P1
-
- Let's open up the matrix multiplication, so we can find the 'X'
- elements as function of the 'W' elements and 'a'.
-
- x0*w0 + x1*w2 = -P0 (1)
- x0*w1 + x1*w3 = 0 (2)
- x2*w0 + x3*w2 = 0 (3)
- x2*w1 + x3*w3 = P1 (4)
-
- Let's bring the X with the smaller index in each equation to the left
- side:
-
- x0 = (-w2/w0)*x1 - P0/w0 (5)
- x0 = (-w3/w1)*x1 (6)
- x2 = (-w2/w0)*x3 (7)
- x2 = (-w3/w1)*x3 + P1/w1 (8)
-
- Using (5) and (6) we can find x0 and x1, by first eliminating x0:
-
- (-w2/w0)*x1 - P0/w0 = (-w3/w1)*x1 -> (w3/w1 - w2/w0) * x1 = P0/w0
-
- For easy reading/writing, let's define: A = (w3/w1 - w2/w0)
-
- --> x1 = P0 / w0 / A
- --> x0 = -1 * x1 * w3 / w1
-
- Similar to the above, we can find x2 and x3 from (7) and (8):
-
- (-w2/w0)*x3 = (-w3/w1)*x3 + P1/w1 -> (w3/w1 - w2/w0) * x3 = P1/w1
-
- --> x3 = P1 / w1 / A
- --> x2 = -1 * x3 * w2 / w0
-
- Note that when the image is already aligned (off-diagonals are zero),
- only the signs of the diagonal elements matter. */
- if( w[1]==0.0f && w[2]==0.0f )
- {
- x[0] = w[0]<0 ? 1.0f : -1.0f; /* Has to be negative. */
- x[1] = 0.0f;
- x[2] = 0.0f;
- x[3] = w[3]>0 ? 1.0f : -1.0f; /* Has to be positive. */
- }
- else if (w[0]==0.0f && w[3]==0.0f )
- {
- x[0] = 0.0f;
- x[1] = w[1]<0 ? 1.0f : -1.0f; /* Has to be negative. */
- x[2] = w[2]>0 ? 1.0f : -1.0f; /* Has to be positive. */
- x[3] = 0.0f;
- }
- else
- {
- A = (w[3]/w[1]) - (w[2]/w[0]);
- x[1] = P[0] / w[0] / A;
- x[3] = P[1] / w[1] / A;
- x[0] = -1 * x[1] * w[3] / w[1];
- x[2] = -1 * x[3] * w[2] / w[0];
-
- /* When the input matrix is flipped, the diagonal elements of the
- necessary rotation will have a different sign. */
- if(x[0]*x[3]<0)
- if(!p->cp.quiet)
- error(0, 0, "%s (hdu %s): WARNING (bug #55217)! The alignment may "
- "not be correct! This is a recognized bug, we just haven't "
- "had enough time to fix it yet, any help or support is "
- "most welcome.\n\n"
- "This may be due to the East and North not being left-hand "
- "oriented. If this is the case, you can fix it by flipping "
- "the image along the problematic axis (with '--flip=1,0' or "
- "'--flip=0,1') and re-running Warp with '--align'. Note "
- "that you can't yet call '--flip' in the same command as "
- "'--align'", p->inputname, p->cp.hdu);
- }
-
- /* For a check:
- printf("ps: (%e, %e)\n", P[0], P[1]);
- printf("w:\n");
- printf(" %.8e %.8e\n", w[0], w[1]);
- printf(" %.8e %.8e\n", w[2], w[3]);
- printf("x:\n");
- printf(" %.8e %.8e\n", x[0], x[1]);
- printf(" %.8e %.8e\n", x[2], x[3]);
- exit(0);
- */
-
- /* Put the matrix elements into the output array: */
- tmatrix[0]=x[0]; tmatrix[1]=x[1]; tmatrix[2]=0.0f;
- tmatrix[3]=x[2]; tmatrix[4]=x[3]; tmatrix[5]=0.0f;
- tmatrix[6]=0.0f; tmatrix[7]=0.0f; tmatrix[8]=1.0f;
-}
-
-
-
-
-
static void
ui_matrix_inplacw_multiply(double *in, double *with)
{
@@ -662,10 +701,6 @@ ui_matrix_from_modular(struct warpparams *p)
structure.*/
switch(pop->status)
{
- case UI_KEY_ALIGN:
- ui_matrix_make_align(p, module);
- break;
-
case UI_KEY_ROTATE:
s = sin( v1 * M_PI / 180 );
c = cos( v1 * M_PI / 180 );
@@ -872,6 +907,10 @@ ui_matrix_finalize(struct warpparams *p)
char *
ui_set_suffix(struct warpparams *p)
{
+ /* Return the suffix as soon as nonlinear mode is detected. Just
+ ignore the rest. */
+ if(p->nonlinearmode) return "_aligned.fits";
+
/* A small independent sanity check: we either need a matrix or at least
one modular warping. */
if(p->matrix==NULL && p->modularll==NULL) ui_error_no_warps();
@@ -881,9 +920,6 @@ ui_set_suffix(struct warpparams *p)
if(p->matrix==NULL && p->modularll->next==NULL)
switch(p->modularll->status)
{
- case UI_KEY_ALIGN:
- return "_aligned.fits";
-
case UI_KEY_ROTATE:
return "_rotated.fits";
@@ -929,8 +965,8 @@ ui_preparations(struct warpparams *p)
p->cp.output=gal_checkset_automatic_output(&p->cp, p->inputname,
ui_set_suffix(p));
- /* Prepare the final warping matrix. */
- ui_matrix_finalize(p);
+ /* Prepare the final warping matrix if in linear mode. */
+ if(p->nonlinearmode == 0) ui_matrix_finalize(p);
}
@@ -959,6 +995,7 @@ void
ui_read_check_inputs_setup(int argc, char *argv[], struct warpparams *p)
{
double *matrix;
+ uint8_t disttype;
struct gal_options_common_params *cp=&p->cp;
@@ -1010,19 +1047,30 @@ ui_read_check_inputs_setup(int argc, char *argv[],
struct warpparams *p)
/* Everything is ready, notify the user of the program starting. */
if(!p->cp.quiet)
{
- matrix=p->matrix->array;
printf(PROGRAM_NAME" "PACKAGE_VERSION" started on %s",
ctime(&p->rawtime));
printf(" Using %zu CPU thread%s\n", p->cp.numthreads,
p->cp.numthreads==1 ? "." : "s.");
printf(" Input: %s (hdu: %s)\n", p->inputname, p->cp.hdu);
- printf(" matrix:"
- "\n\t%.4f %.4f %.4f"
- "\n\t%.4f %.4f %.4f"
- "\n\t%.4f %.4f %.4f\n",
- matrix[0], matrix[1], matrix[2],
- matrix[3], matrix[4], matrix[5],
- matrix[6], matrix[7], matrix[8]);
+
+ if(p->nonlinearmode)
+ {
+ disttype=gal_wcs_distortion_identify(p->input->wcs);
+ if(disttype!=GAL_WCS_DISTORTION_INVALID)
+ printf(" matrix: '%s' distortion from WCS of input.\n",
+ gal_wcs_distortion_to_string(disttype));
+ }
+ else
+ {
+ matrix=p->matrix->array;
+ printf(" matrix:"
+ "\n\t% .4f % .4f % .4f"
+ "\n\t% .4f % .4f % .4f"
+ "\n\t% .4f % .4f % .4f\n",
+ matrix[0], matrix[1], matrix[2],
+ matrix[3], matrix[4], matrix[5],
+ matrix[6], matrix[7], matrix[8]);
+ }
}
}
@@ -1052,12 +1100,23 @@ void
ui_free_report(struct warpparams *p, struct timeval *t1)
{
/* Free the allocated arrays: */
+ if(p->inverse) free(p->inverse);
+ if(p->wa.cdelt) p->wa.cdelt=NULL;
+ if(p->cdelt) gal_data_free(p->cdelt);
+ if(p->matrix) gal_data_free(p->matrix);
+ if(p->inwcsmatrix) free(p->inwcsmatrix);
+ if(p->modularll) gal_data_free(p->modularll);
+
+ if(p->wa.input) p->wa.input=NULL;
+ if(p->wa.ctype) gal_data_free(p->wa.ctype);
+ if(p->wa.center) gal_data_free(p->wa.center);
+ if(p->wa.widthinpix) gal_data_free(p->wa.widthinpix);
+
free(p->cp.hdu);
free(p->cp.output);
gal_data_free(p->input);
- gal_data_free(p->matrix);
- if(p->pixelscale) free(p->pixelscale);
- if(p->inwcsmatrix) free(p->inwcsmatrix);
+ gal_data_free(p->output);
+
/* Report how long the operation took. */
if(!p->cp.quiet)
diff --git a/bin/warp/ui.h b/bin/warp/ui.h
index 763e613c..f976e6cb 100644
--- a/bin/warp/ui.h
+++ b/bin/warp/ui.h
@@ -5,6 +5,7 @@ Warp is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <mohammad@akhlaghi.org>
Contributing author(s):
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
Copyright (C) 2016-2022 Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
@@ -33,7 +34,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Option groups particular to this program. */
enum program_args_groups
{
- UI_GROUP_WARPS = GAL_OPTIONS_GROUP_AFTER_COMMON,
+ UI_GROUP_ALIGN = GAL_OPTIONS_GROUP_AFTER_COMMON,
+ UI_GROUP_WARPS
};
@@ -42,7 +44,7 @@ enum program_args_groups
/* Available letters for short options:
- b d g i j l n u v w x y z
+ a b d g i j l n u v z
A B E G H J L O Q R W X Y
*/
enum option_keys_enum
@@ -50,7 +52,6 @@ enum option_keys_enum
/* With short-option version. */
UI_KEY_KEEPWCS = 'k',
UI_KEY_COVEREDFRAC = 'C',
- UI_KEY_ALIGN = 'a',
UI_KEY_ROTATE = 'r',
UI_KEY_SCALE = 's',
UI_KEY_FLIP = 'f',
@@ -58,12 +59,18 @@ enum option_keys_enum
UI_KEY_TRANSLATE = 't',
UI_KEY_PROJECT = 'p',
UI_KEY_MATRIX = 'm',
- UI_KEY_CENTERONCORNER = 'c',
+ UI_KEY_CDELT = 'x',
+ UI_KEY_INTERPSAMPLING = 'y',
+ UI_KEY_CENTER = 'c',
+ UI_KEY_WIDTHINPIX = 'w',
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
- UI_KEY_HSTARTWCS = 1000,
+ UI_KEY_CENTERONCORNER = 1000,
+ UI_KEY_EDGESAMPLING,
+ UI_KEY_HSTARTWCS,
UI_KEY_HENDWCS,
+ UI_KEY_CTYPE,
};
diff --git a/bin/warp/warp.c b/bin/warp/warp.c
index e328efb0..1c0ad5a2 100644
--- a/bin/warp/warp.c
+++ b/bin/warp/warp.c
@@ -2,9 +2,10 @@
Warp - Warp images using projective mapping.
Warp is part of GNU Astronomy Utilities (Gnuastro) package.
-Original author:
+Corresponding author:
Mohammad Akhlaghi <mohammad@akhlaghi.org>
Contributing author(s):
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
Copyright (C) 2015-2022 Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
@@ -30,8 +31,12 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+#include <gnuastro/warp.h>
#include <gnuastro/polygon.h>
#include <gnuastro/pointer.h>
+#include <gnuastro/threads.h>
+
+#include <gnuastro-internal/timing.h>
#include "main.h"
#include "warp.h"
@@ -40,15 +45,13 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-
-
/***************************************************************/
/************** MACROS ******************/
/***************************************************************/
-/* Multiply a 2 element vector with a transformation matrix and put
- the result in the 2 element output array. It is assumed that the
- input is from a flat coordinate systemy. */
-#define mappoint(V, T, O) \
+/* Multiply a 2 element vector with a transformation matrix and put the
+ result in the 2 element output array. It is assumed that the input is
+ from a flat coordinate systemy. */
+#define WARP_MAPPOINT(V, T, O) \
{ \
(O)[0]=( ( (T)[0]*(V)[0] + (T)[1]*(V)[1] + (T)[2] ) \
/ ( (T)[6]*(V)[0] + (T)[7]*(V)[1] + (T)[8] ) ); \
@@ -60,39 +63,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-/* A pixel's center is an integer value. This function will give the
- integer value that is nearest to a floating point number. Works for
- both positive and negative values and includes floating point
- errors.
-
- nearestint_halfhigher(0.5f) --> 1.0f
-*/
-#define nearestint_halfhigher(D) \
- (ceil((D)) - (D) > 0.5f \
- + GAL_POLYGON_ROUND_ERR ? ceil((D))-1.0f : ceil((D)))
-
-
-
-
-
-/* Similar to 'nearestint_halfhigher' but:
-
- nearestint_halflower(0.5f) --> 0.0f;
- */
-#define nearestint_halflower(D) \
- (ceil((D)) - (D) > 0.5f - GAL_POLYGON_ROUND_ERR ? ceil((D))-1.0f : ceil((D)))
-
-
-
-
-
-#define ceilwitherr(D) \
- ( fabs( nearbyint((D)) - (D) ) < GAL_POLYGON_ROUND_ERR \
- ? nearbyint((D)) : nearbyint((D))+1 )
-
-
-
-
@@ -113,10 +83,10 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/************** Processing function ******************/
/***************************************************************/
static void *
-warp_onthread(void *inparam)
+warp_onthread_linear(void *inparam)
{
- struct iwpparams *iwp=(struct iwpparams*)inparam;
- struct warpparams *p=iwp->p;
+ struct gal_threads_params *tprm=(struct gal_threads_params *)inparam;
+ struct warpparams *p=(struct warpparams *)tprm->params;
size_t *extinds=p->extinds, *ordinds=p->ordinds;
long is0=p->input->dsize[0], is1=p->input->dsize[1];
@@ -126,8 +96,10 @@ warp_onthread(void *inparam)
double ocrn[8], icrn_base[8], icrn[8], *output=p->output->array;
double pcrn[8], *outfpixval=p->outfpixval, ccrn[GAL_POLYGON_MAX_CORNERS];
- for(i=0; (ind=iwp->indexs[i])!=GAL_BLANK_SIZE_T; ++i)
+ for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
{
+ ind=tprm->indexs[i];
+
/* Initialize the output pixel value: */
numinput=0;
output[ind]=filledarea=0.0f;
@@ -148,14 +120,19 @@ warp_onthread(void *inparam)
/* Transform the four corners of the output pixel to the input
image coordinates. */
- for(j=0;j<4;++j) mappoint(&ocrn[j*2], p->inverse, &icrn_base[j*2]);
+ for(j=0;j<4;++j) WARP_MAPPOINT(&ocrn[j*2], p->inverse,
+ &icrn_base[j*2]);
/* Using the known relationships between the vertice locations,
put everything in the right place: */
- xstart = nearestint_halfhigher( icrn_base[extinds[0]] );
- xend = nearestint_halflower( icrn_base[extinds[1]] ) + 1;
- ystart = nearestint_halfhigher( icrn_base[extinds[2]] );
- yend = nearestint_halflower( icrn_base[extinds[3]] ) + 1;
+ xstart =
+ GAL_DIMENSION_NEARESTINT_HALFHIGHER( icrn_base[extinds[0]] );
+ xend =
+ GAL_DIMENSION_NEARESTINT_HALFLOWER( icrn_base[extinds[1]] ) + 1;
+ ystart =
+ GAL_DIMENSION_NEARESTINT_HALFHIGHER( icrn_base[extinds[2]] );
+ yend =
+ GAL_DIMENSION_NEARESTINT_HALFLOWER( icrn_base[extinds[3]] ) + 1;
icrn[0]=icrn_base[ordinds[0]*2]; icrn[1]=icrn_base[ordinds[0]*2+1];
icrn[2]=icrn_base[ordinds[1]*2]; icrn[3]=icrn_base[ordinds[1]*2+1];
icrn[4]=icrn_base[ordinds[2]*2]; icrn[5]=icrn_base[ordinds[2]*2+1];
@@ -250,12 +227,8 @@ warp_onthread(void *inparam)
if(numinput==0) output[ind]=NAN;
}
-
-
- /* Wait until all other threads finish. */
- if(p->cp.numthreads>1)
- pthread_barrier_wait(iwp->b);
-
+ /* Wait for all the other threads to finish, then return. */
+ if(tprm->b) { pthread_barrier_wait(tprm->b); }
return NULL;
}
@@ -275,7 +248,6 @@ warp_onthread(void *inparam)
-
/***************************************************************/
/************** Preparations ******************/
/***************************************************************/
@@ -293,7 +265,7 @@ warp_onthread(void *inparam)
the image altough the scale might change.
*/
static void
-warp_preparations(struct warpparams *p)
+warp_linear_init(struct warpparams *p)
{
double is0=p->input->dsize[0], is1=p->input->dsize[1];
@@ -306,18 +278,20 @@ warp_preparations(struct warpparams *p)
double input[8]={ 0.5f, 0.5f, is1+0.5f, 0.5f,
0.5f, is0+0.5f, is1+0.5f, is0+0.5f };
+
/* Find the range of pixels of the input image. All the input positions
are moved to the negative by half a pixel since the center of the
pixel is an integer value.*/
for(i=0;i<4;++i)
{
- mappoint(&input[i*2], matrix, &output[i*2]);
+ WARP_MAPPOINT(&input[i*2], matrix, &output[i*2]);
if(output[i*2]<xmin) xmin = output[i*2];
if(output[i*2]>xmax) xmax = output[i*2];
if(output[i*2+1]<ymin) ymin = output[i*2+1];
if(output[i*2+1]>ymax) ymax = output[i*2+1];
}
+
/* For a check:
for(i=0;i<4;++i)
printf("(%.3f, %.3f) --> (%.3f, %.3f)\n",
@@ -327,14 +301,18 @@ warp_preparations(struct warpparams *p)
xmin, xmax, ymin, ymax);
*/
+
/* Set the final size of the image. The X axis is horizontal. The reason
we are using the halflower variation of 'nearestint' for the maximums
is that these points are the farthest extremes of the input image. If
they are half a pixel value, they should point to the pixel before. */
- dsize[1]=nearestint_halflower(xmax)-nearestint_halfhigher(xmin)+1;
- dsize[0]=nearestint_halflower(ymax)-nearestint_halfhigher(ymin)+1;
- p->outfpixval[0]=nearestint_halfhigher(xmin);
- p->outfpixval[1]=nearestint_halfhigher(ymin);
+ dsize[1] = GAL_DIMENSION_NEARESTINT_HALFLOWER(xmax) \
+ - GAL_DIMENSION_NEARESTINT_HALFHIGHER(xmin)+1;
+ dsize[0] = GAL_DIMENSION_NEARESTINT_HALFLOWER(ymax) \
+ - GAL_DIMENSION_NEARESTINT_HALFHIGHER(ymin)+1;
+ p->outfpixval[0]=GAL_DIMENSION_NEARESTINT_HALFHIGHER(xmin);
+ p->outfpixval[1]=GAL_DIMENSION_NEARESTINT_HALFHIGHER(ymin);
+
/* If we have translation, the 'dsize's and 'outfpixval's should be
corrected. Note that centeroncorner is also a translation operation,
@@ -347,6 +325,7 @@ warp_preparations(struct warpparams *p)
if(ymin>0) p->outfpixval[1]=0;
}
+
/* For a check:
printf("Wrapped:\n");
printf("dsize [C]: (%zu, %zu)\n", dsize[0], dsize[1]);
@@ -354,6 +333,7 @@ warp_preparations(struct warpparams *p)
p->outfpixval[1]);
*/
+
/* We now know the size of the output and the starting and ending
coordinates in the output image (bottom left corners of pixels)
for the transformation. */
@@ -371,7 +351,7 @@ warp_preparations(struct warpparams *p)
{
ocrn[i*2] += p->outfpixval[0];
ocrn[i*2+1] += p->outfpixval[1];
- mappoint(&ocrn[i*2], p->inverse, &icrn[i*2]);
+ WARP_MAPPOINT(&ocrn[i*2], p->inverse, &icrn[i*2]);
}
@@ -403,6 +383,7 @@ warp_preparations(struct warpparams *p)
if(icrn[i*2+1]>ymax) { ymax=icrn[i*2+1]; extinds[3]=i*2+1; }
}
+
/* For a check:
for(i=0;i<4;++i)
printf("(%.3f, %.3f) --> (%.3f, %.3f)\n",
@@ -416,18 +397,51 @@ warp_preparations(struct warpparams *p)
+static void
+warp_write_to_file(struct warpparams *p, int hasmatrix)
+{
+ size_t i;
+ char keyword[9*FLEN_KEYWORD];
+ gal_fits_list_key_t *headers=NULL;
+ double *m= hasmatrix ? p->matrix->array : NULL;
+
+ /* Add the appropriate headers: */
+ gal_fits_key_write_filename("INF", p->inputname, &headers,
+ 0, p->cp.quiet);
+ if(hasmatrix)
+ for(i=0;i<9;++i)
+ {
+ sprintf(&keyword[i*FLEN_KEYWORD], "WMTX%zu_%zu", i/3+1, i%3+1);
+ gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64,
+ &keyword[i*FLEN_KEYWORD], 0, &m[i], 0,
+ "Warp matrix element value", 0, NULL, 0);
+ }
+
+ /* Save the output into the proper type and write it. */
+ if(p->cp.type && p->cp.type!=p->output->type)
+ p->output=gal_data_copy_to_new_type_free(p->output, p->cp.type);
+ gal_fits_img_write(p->output, p->cp.output, headers, PROGRAM_NAME);
+
+ /* Write the configuration keywords. */
+ gal_fits_key_write_filename("input", p->inputname, &p->cp.okeys,
+ 1, p->cp.quiet);
+ gal_fits_key_write_config(&p->cp.okeys, "Warp configuration",
+ "WARP-CONFIG", p->cp.output, "0");
+}
+
+
+
+
+
/* Correct the WCS coordinates (Multiply the 2x2 PC matrix of the WCS
structure by the INVERSE of the transform in 2x2). Then Multiply the
crpix array with the ACTUAL transformation matrix. */
void
-correct_wcs_save_output(struct warpparams *p)
+warp_write_wcs_linear(struct warpparams *p)
{
- size_t i;
- double tcrpix[3];
- char keyword[9*FLEN_KEYWORD];
double *m=p->matrix->array, diff;
struct wcsprm *wcs=p->output->wcs;
- gal_fits_list_key_t *headers=NULL;
+ double tcrpix[3], *ps=p->cdelt->array;
double *crpix=wcs?wcs->crpix:NULL, *w=p->inwcsmatrix;
/* 'tinv' is the 2 by 2 inverse matrix. Recall that 'p->inverse' is 3 by
@@ -469,32 +483,19 @@ correct_wcs_save_output(struct warpparams *p)
if( fabs(wcs->pc[1])<ABSOLUTEFLTERROR ) wcs->pc[1]=0.0f;
if( fabs(wcs->pc[2])<ABSOLUTEFLTERROR ) wcs->pc[2]=0.0f;
diff=fabs(wcs->pc[0])-fabs(wcs->pc[3]);
- if( fabs(diff/p->pixelscale[0])<RELATIVEFLTERROR )
- wcs->pc[3]=( (wcs->pc[3] < 0.0f ? -1.0f : 1.0f) * fabs(wcs->pc[0]) );
+ if( fabs(diff/ps[0])<RELATIVEFLTERROR )
+ wcs->pc[3]=( (wcs->pc[3] < 0.0f ? -1.0f : 1.0f)
+ * fabs(wcs->pc[0]) );
}
- /* Add the appropriate headers: */
- gal_fits_key_write_filename("INF", p->inputname, &headers,
- 0, p->cp.quiet);
- for(i=0;i<9;++i)
- {
- sprintf(&keyword[i*FLEN_KEYWORD], "WMTX%zu_%zu", i/3+1, i%3+1);
- gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64,
- &keyword[i*FLEN_KEYWORD], 0, &m[i], 0,
- "Warp matrix element value", 0, NULL, 0);
- }
+ /* Write the final keywords and the file. */
+ warp_write_to_file(p, 1);
+}
+
+
+
- /* Save the output into the proper type and write it. */
- if(p->cp.type && p->cp.type!=p->output->type)
- p->output=gal_data_copy_to_new_type_free(p->output, p->cp.type);
- gal_fits_img_write(p->output, p->cp.output, headers, PROGRAM_NAME);
- /* Write the configuration keywords. */
- gal_fits_key_write_filename("input", p->inputname, &p->cp.okeys,
- 1, p->cp.quiet);
- gal_fits_key_write_config(&p->cp.okeys, "Warp configuration",
- "WARP-CONFIG", p->cp.output, "0");
-}
@@ -521,80 +522,51 @@ correct_wcs_save_output(struct warpparams *p)
void
warp(struct warpparams *p)
{
- int err;
- pthread_t t; /* All thread ids saved in this, not used. */
- char *mmapname;
- pthread_attr_t attr;
- pthread_barrier_t b;
- struct iwpparams *iwp;
- size_t nt=p->cp.numthreads;
- size_t i, nb, *indexs, thrdcols;
-
-
- /* Array keeping thread parameters for each thread. */
- errno=0;
- iwp=malloc(nt*sizeof *iwp);
- if(iwp==NULL)
- error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for iwp",
- __func__, nt*sizeof *iwp);
-
+ struct timeval t0;
+ gal_warp_wcsalign_t *wa=&p->wa;
- /* Prepare the output array and all the necessary things: */
- warp_preparations(p);
-
-
- /* Distribute the output pixels into the threads: */
- mmapname=gal_threads_dist_in_threads(p->output->size, nt,
- p->cp.minmapsize, p->cp.quietmmap,
- &indexs, &thrdcols);
-
-
- /* Start the warp. */
- if(nt==1)
+ /* Do the preparations and set the pointers to the functions to use. */
+ if( p->nonlinearmode )
{
- iwp[0].p=p;
- iwp[0].indexs=indexs;
- warp_onthread(&iwp[0]);
+ /* Calculate and allocate the output image size and WCS */
+ if(!p->cp.quiet)
+ {
+ gal_timing_report(NULL, "Initializing the output image...", 1);
+ gettimeofday(&t0, NULL);
+ }
+ gal_warp_wcsalign_init(wa);
+
+ /* Fill the output image */
+ if(!p->cp.quiet)
+ {
+ gal_timing_report(&t0, "Done", 2);
+ gal_timing_report(NULL, "Warping the input image...", 1);
+ gettimeofday(&t0, NULL);
+ }
+ gal_threads_spin_off(gal_warp_wcsalign_onthread, wa,
+ wa->output->size, wa->numthreads,
+ wa->input->minmapsize, wa->input->quietmmap);
+ if(!p->cp.quiet) gal_timing_report(&t0, "Done", 2);
+ p->output=wa->output;
+ wa->output=NULL; /* must be here! */
+ gal_warp_wcsalign_free(wa);
+
+ /* Write the final keywords and the file. */
+ warp_write_to_file(p, 0);
}
else
{
- /* Initialize the attributes. Note that this running thread
- (that spinns off the nt threads) is also a thread, so the
- number the barrier should be one more than the number of
- threads spinned off. */
- if(p->output->size<nt) nb=p->output->size+1;
- else nb=nt+1;
- gal_threads_attr_barrier_init(&attr, &b, nb);
-
- /* Spin off the threads: */
- for(i=0;i<nt;++i)
- if(indexs[i*thrdcols]!=GAL_BLANK_SIZE_T)
- {
- iwp[i].p=p;
- iwp[i].b=&b;
- iwp[i].indexs=&indexs[i*thrdcols];
- err=pthread_create(&t, &attr, warp_onthread, &iwp[i]);
- if(err)
- error(EXIT_FAILURE, 0, "%s: can't create thread %zu",
- __func__, i);
- }
-
- /* Wait for all threads to finish and free the spaces. */
- pthread_barrier_wait(&b);
- pthread_attr_destroy(&attr);
- pthread_barrier_destroy(&b);
- }
+ warp_linear_init(p);
+ /* Fill the output image */
+ gal_threads_spin_off(warp_onthread_linear, p, p->output->size,
+ p->cp.numthreads, p->cp.minmapsize,
+ p->cp.quietmmap);
- /* Save the output. */
- correct_wcs_save_output(p);
- if(!p->cp.quiet)
- printf(" Output: %s\n", p->cp.output);
+ /* Fix the linear matrix before saving the output image to disk */
+ warp_write_wcs_linear(p);
+ }
- /* Free the allocated spaces, note that 'indexs' may be memory-mapped. */
- free(iwp);
- gal_data_free(p->output);
- if(mmapname) gal_pointer_mmap_free(&mmapname, p->cp.quietmmap);
- else free(indexs);
+ if(!p->cp.quiet) { printf(" Output: %s\n", p->cp.output); }
}
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0674db16..467711c6 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -552,11 +552,16 @@ Frequency domain and Fourier operations
Warp
-* Warping basics:: Basics of coordinate transformation.
+* Linear warping basics:: Basics of coordinate transformation.
* Merging multiple warpings:: How to merge multiple matrices.
* Resampling:: Warping an image is re-sampling it.
* Invoking astwarp:: Arguments and options for Warp.
+Invoking Warp
+
+* Align pixels with WCS and remove distortions:: Default operation.
+* Linear warps to be called explicitly:: Other waprs.
+
Data analysis
* Statistics:: Calculate dataset statistics.
@@ -570,7 +575,7 @@ Statistics
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
-* Least squares fitting::
+* Least squares fitting:: Fitting with various parametric functions.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@@ -786,6 +791,7 @@ Gnuastro library
* Labeled datasets:: Working with Segmented/labeled datasets.
* Convolution functions:: Library functions to do convolution.
* Interpolation:: Interpolate (over blank values possibly).
+* Warp library:: Warp pixel grid to a new one.
* Color functions:: Definitions and operations related to colors.
* Git wrappers:: Wrappers for functions in libgit2.
* Python interface:: Functions to help in writing Python wrappers.
@@ -835,11 +841,6 @@ File input output
* EPS files:: Writing to EPS files.
* PDF files:: Writing to PDF files.
-World Coordinate System (@file{wcs.h})
-
-* Freeing a WCS structure:: Demo on creating and properly freeing.
-* WCS Function and global constants:: All functions in gnuastro/wcs.h.
-
Tessellation library (@file{tile.h})
* Independent tiles:: Work on or check independent tiles.
@@ -851,6 +852,7 @@ Library demo programs
* Library demo - inspecting neighbors:: Inspect the neighbors of a pixel.
* Library demo - multi-threaded operation:: Doing an operation on threads.
* Library demo - reading and writing table columns:: Simple Column I/O.
+* Library demo - resampling an image to a new grid::
Developing
@@ -12438,7 +12440,7 @@ However, the WCS coordinate system of
@file{img-eq.fits} will be the equatorial
Fits will automatically extract the current coordinate system of
@file{image.fits} and as long as it is one of the recognized coordinate systems
listed below, it will do the conversion.
@example
-$ astfits image.fits --coordsys=eq-j2000 --output=img-eq.fits
+$ astfits image.fits --wcscoordsys=eq-j2000 --output=img-eq.fits
@end example
The currently recognized coordinate systems are listed below (the most common
one today is @code{eq-j2000}):
@@ -19188,7 +19190,7 @@ Here we specifically mean the pixel grid transformation
which is better conveyed
@cindex Gravitational lensing
Image wrapping is a very important step in astronomy, both in observational
data analysis and in simulating modeled images.
-In modeling, warping an image is necessary when we want to apply grid
transformations to the initial models, for example in simulating gravitational
lensing (Radial warpings are not yet included in Warp).
+In modeling, warping an image is necessary when we want to apply grid
transformations to the initial models, for example in simulating gravitational
lensing.
Observational reasons for warping an image are listed below:
@itemize
@@ -19219,7 +19221,7 @@ To do that, we need all the images to be on the same
pixel grid.
@cindex Optical distortion
@cindex Distortion, optical
@item
-@strong{Optical distortion:} (Not yet included in Warp) In wide field images,
the optical distortion that occurs on the outer parts of the focal plane will
make accurate comparison of the objects at various locations impossible.
+@strong{Optical distortion:} In wide field images, the optical distortion that
occurs on the outer parts of the focal plane will make accurate comparison of
the objects at various locations impossible.
It is therefore necessary to warp the image and correct for those distortions
prior to the analysis.
@cindex ACS
@@ -19235,14 +19237,14 @@ It is therefore necessary to warp the image and
correct for those distortions pr
@end itemize
@menu
-* Warping basics:: Basics of coordinate transformation.
+* Linear warping basics:: Basics of coordinate transformation.
* Merging multiple warpings:: How to merge multiple matrices.
* Resampling:: Warping an image is re-sampling it.
* Invoking astwarp:: Arguments and options for Warp.
@end menu
-@node Warping basics, Merging multiple warpings, Warp, Warp
-@subsection Warping basics
+@node Linear warping basics, Merging multiple warpings, Warp, Warp
+@subsection Linear warping basics
@cindex Scaling
@cindex Coordinate transformation
@@ -19353,7 +19355,7 @@ A short but instructive and illustrated review of
affine, projective and also bi
Paul S. Heckbert. 1989. @emph{Fundamentals of Texture mapping and Image
Warping}, Master's thesis at University of California, Berkeley.
Note that since points are defined as row vectors there, the matrix is the
transpose of the one discussed here.}.
-@node Merging multiple warpings, Resampling, Warping basics, Warp
+@node Merging multiple warpings, Resampling, Linear warping basics, Warp
@subsection Merging multiple warpings
@cindex Commutative property
@@ -19361,7 +19363,7 @@ Note that since points are defined as row vectors
there, the matrix is the trans
@cindex Multiplication, matrix
@cindex Non-commutative operations
@cindex Operations, non-commutative
-In @ref{Warping basics} we saw how a basic warp/transformation can be
represented with a matrix.
+In @ref{Linear warping basics} we saw how a basic warp/transformation can be
represented with a matrix.
To make more complex warpings (for example to define a translation, rotation
and scale as one warp) the individual matrices have to be multiplied through
matrix multiplication.
However matrix multiplication is not commutative, so the order of the set of
matrices you use for the multiplication is going to be very important.
@@ -19396,22 +19398,26 @@ These three operations can be merged in one operation
by calculating the matrix
A digital image is composed of discrete `picture elements' or `pixels'.
When a real image is created from a camera or detector, each pixel's area is
used to store the number of photo-electrons that were created when incident
photons collided with that pixel's surface area.
This process is called the `sampling' of a continuous or analog data into
digital data.
-When we change the pixel grid of an image or warp it as we defined in
@ref{Warping basics}, we have to `guess' the flux value of each pixel on the
new grid based on the old grid, or re-sample it.
-Because of the `guessing', any form of warping on the data is going to degrade
the image and mix the original pixel values with each other.
+
+When we change the pixel grid of an image, or ``warp'' it, we have to
calculate the flux value of each pixel on the new grid based on the old grid,
or re-sample it.
+Because of the calculation (as opposed to observation), any form of warping on
the data is going to degrade the image and mix the original pixel values with
each other.
So if an analysis can be done on an unwarped data image, it is best to leave
the image untouched and pursue the analysis.
-However as discussed in @ref{Warp} this is not possible most of the times, so
we have to accept the problem and re-sample the image.
+However as discussed in @ref{Warp} this is not possible in some scenarios and
re-sampling is necessary.
@cindex Point pixels
@cindex Interpolation
+@cindex Sampling theorem
@cindex Bicubic interpolation
@cindex Signal to noise ratio
@cindex Bi-linear interpolation
@cindex Interpolation, bicubic
@cindex Interpolation, bi-linear
-In most applications of image processing, it is sufficient to consider each
pixel to be a point and not an area.
-This assumption can significantly speed up the processing of an image and also
the simplicity of the code.
-It is a fine assumption when the signal to noise ratio of the objects are very
large.
-The question will then be one of interpolation because you have multiple
points distributed over the output image and you want to find the values at the
pixel centers.
+When the FWHM of the PSF of the camera is much larger than the pixel scale
(see @ref{Sampling theorem}) we are sampling the signal in a much higher
resolution than the camera can offer.
+This is usually the case in many applications of image processing
(non-astronomical imaging).
+In such cases, we can consider each pixel to be a point and not an area: the
PSF doesn't vary much over a single pixel.
+
+Approximating a pixel's area to a point can significantly speed up the
resampling and also the simplicity of the code.
+Because resampling becomes a problem of interpolation: points of the input
grid need to be interpolated at certain other points (over the output grid).
To increase the accuracy, you might also sample more than one point from
within a pixel giving you more points for a more accurate interpolation in the
output grid.
@cindex Image edges
@@ -19420,12 +19426,12 @@ However, interpolation has several problems.
The first one is that it will depend on the type of function you want to
assume for the interpolation.
For example you can choose a bi-linear or bi-cubic (the `bi's are for the 2
dimensional nature of the data) interpolation method.
For the latter there are various ways to set the constants@footnote{see
@url{http://entropymine.com/imageworsener/bicubic/} for a nice introduction.}.
-Such functional interpolation functions can fail seriously on the edges of an
image.
-They will also need normalization so that the flux of the objects before and
after the warpings are comparable.
-The most basic problem with such techniques is that they are based on a point
while a detector pixel is an area.
-They add a level of subjectivity to the data (make more assumptions through
the functions than the data can handle).
-For most applications this is fine, but in scientific applications where
detection of the faintest possible galaxies or fainter parts of bright galaxies
is our aim, we cannot afford this loss.
-Because of these reasons Warp will not use such interpolation techniques.
+Such parametric interpolation functions can fail seriously on the edges of an
image, or when there is a sharp change in value (for example the bleeding
saturation of bright stars in astronomical CCDs).
+They will also need normalization so that the flux of the objects before and
after the warping is comparable.
+
+The parametric nature of these methods adds a level of subjectivity to the
data (it makes more assumptions through the functions than the data can handle).
+For most applications this is fine (as discussed above: when the PSF is
over-sampled), but in scientific applications where we push our instruments to
the limit and the aim is the detection of the faintest possible galaxies or
fainter parts of bright galaxies, we cannot afford this loss.
+Because of these reasons Warp will not use parametric interpolation techniques.
@cindex Drizzle
@cindex Pixel mixing
@@ -19434,20 +19440,33 @@ Warp will do interpolation based on ``pixel
mixing''@footnote{For a graphic demo
This is also what the Hubble Space Telescope pipeline calls
``Drizzling''@footnote{@url{http://en.wikipedia.org/wiki/Drizzle_(image_processing)}}.
This technique requires no functions, it is thus non-parametric.
It is also the closest we can get (make least assumptions) to what actually
happens on the detector pixels.
-The basic idea is that you reverse-transform each output pixel to find which
pixels of the input image it covers and what fraction of the area of the input
pixels are covered.
-To find the output pixel value, you simply sum the value of each input pixel
weighted by the overlapfraction (between 0 to 1) of the output pixel and that
input pixel.
-Through this process, pixels are treated as an area not as a point (which is
how detectors create the image), also the brightness (see @ref{Brightness flux
magnitude}) of an object will be left completely unchanged.
-If there are very high spatial-frequency signals in the image (for example
fringes) which vary on a scale smaller than your output image pixel size, pixel
mixing can cause
ailiasing@footnote{@url{http://en.wikipedia.org/wiki/Aliasing}}.
-So if the input image has fringes, they have to be calculated and removed
separately (which would naturally be done in any astronomical application).
-Because of the PSF no astronomical target has a sharpchange in the signal so
this issue is less important for astronomical applications, see @ref{PSF}.
+In pixel mixing, the basic idea is that you reverse-transform each output
pixel to find which pixels of the input image it covers, and what fraction of
the area of the input pixels are covered by that output pixel.
+We then multiply each input pixel's value by the fraction of its area that
overlaps with the output pixel (between 0 to 1).
+The output's pixel value is derived by summing all these multiplications for
the input pixels that it covers.
+
+Through this process, pixels are treated as an area not as a point (which is
how detectors create the image), also the brightness (see @ref{Brightness flux
magnitude}) of an object will be fully preserved.
+Since it involves the mixing of the input's pixel values, this pixel mixing
method is a form of @ref{Spatial domain convolution}.
+Therefore, after comparing the input and output, you will notice that the
output is slightly smoothed, thus boosting the more diffuse signal, but
creating correlated noise.
+In astronomical imaging the correlated noise will be decreased later when you
stack many exposures.
+
+If there are very high spatial-frequency signals in the image (for example
fringes) which vary on a scale @emph{smaller than} your output image pixel size
(this is rarely the case in astronomical imaging), pixel mixing can cause
ailiasing@footnote{@url{http://en.wikipedia.org/wiki/Aliasing}}.
+Therefore, in case such fringes are present, they have to be calculated and
removed separately (which would naturally be done in any astronomical reduction
pipeline).
+Because of the PSF, no astronomical target has a sharp change in their signal.
+Thus this issue is less important for astronomical applications, see @ref{PSF}.
+To find the overlap area of the output pixel over the input pixels, we need to
define polygons and clip them (find the overlap).
+Usually, it is sufficient to define a pixel with a four-vertice polygon.
+However, when a non-linear distortion (for example @code{SIP} or @code{TPV})
is present and the distortion is significant over an output pixel's size
(usually far from the reference point), the shadow of the output pixel on the
input grid can be curved.
+To account for such cases (which can only happen when correcting for
non-linear distortions), Warp has the @option{--edgesampling} option to sample
the output pixel over more vertices.
+For more, see the description of this option in @ref{Align pixels with WCS and
remove distortions}.
@node Invoking astwarp, , Resampling, Warp
@subsection Invoking Warp
-Warp an input dataset into a new grid.
-Any homographic warp (for example scaling, rotation, translation, projection)
is acceptable, see @ref{Warping basics} for the definitions.
+Warp an input image into a new pixel grid by pixel mixing (see
@ref{Resampling}).
+Without any options, Warp will remove any non-linear distortions from the
image and align the output pixel coordinates to its WCS coordinates.
+Any homographic warp (for example scaling, rotation, translation, projection,
see @ref{Linear warping basics}) can also be done by calling the relevant
option explicitly.
The general template for invoking Warp is:
@example
@@ -19458,15 +19477,26 @@ $ astwarp [OPTIONS...] InputImage
One line examples:
@example
+## Align image with celestial coordinates and remove any distortion
+$ astwarp image.fits
+
+## Align four exposures to same pixel grid and stack them with
+## Arithmetic program's sigma-clipped mean operator (out of many
+## stacking operators, see Arithmetic's documentation).
+$ grid="--center=1.234,5.678 --widthinpix=1001,1001 --cdelt=0.2/3600"
+$ astwarp a.fits $grid --output=A.fits
+$ astwarp b.fits $grid --output=B.fits
+$ astwarp c.fits $grid --output=C.fits
+$ astwarp d.fits $grid --output=D.fits
+$ astarithmetic A.fits B.fits C.fits D.fits 4 5 0.2 sigclip-mean \
+ -g1 --output=stack.fits
+
## Rotate and then scale input image:
$ astwarp --rotate=37.92 --scale=0.8 image.fits
## Scale, then translate the input image:
$ astwarp --scale 8/3 --translate 2.1 image.fits
-## Align raw image with celestial coordinates:
-$ astwarp --align rawimage.fits --output=aligned.fits
-
## Directly input a custom warping matrix (using fraction):
$ astwarp --matrix=1/5,0,4/10,0,1/5,4/10,0,0,1 image.fits
@@ -19474,14 +19504,242 @@ $ astwarp --matrix=1/5,0,4/10,0,1/5,4/10,0,0,1
image.fits
$ astwarp --matrix="0.7071,-0.7071, 0.7071,0.7071" image.fits
@end example
-If any processing is to be done, Warp can accept one file as input.
+If any processing is to be done, Warp needs to be given a 2D FITS image.
As in all Gnuastro programs, when an output is not explicitly set with the
@option{--output} option, the output filename will be set automatically based
on the operation, see @ref{Automatic output}.
For the full list of general options to all Gnuastro programs (including
Warp), please see @ref{Common options}.
-To be the most accurate, the input image will be read as a 64-bit double
precision floating point dataset and all internal processing is done in this
format (including the raw output type).
-You can use the common @option{--type} option to write the output in any type
you want, see @ref{Numeric data types}.
+Warp uses pixel mixing to derive the pixel values of the output image, see
@ref{Resampling}.
+To be the most accurate, the input image will be read as a 64-bit double
precision floating point dataset and all internal processing is done in this
format.
+Upon writing, by default it will be converted to 32-bit single precision
floating point type (actual observational data rarely have such precision!).
+In case you want a different output type, you can use the @option{--type}
option that is common to several Gnuastro programs.
+For example, if your input is a mock image without noise, and you want to
preserve the 64-bit precision, use (with @option{--type=float64}.
+Just note that the file size will also be double!
+For more on the precision of various types, see @ref{Numeric data types}.
+
+By default (if no linear operation is requested), Warp will align the pixel
grid of the input image to the WCS coordinates it contains.
+This operation and the the options that govern it are described in @ref{Align
pixels with WCS and remove distortions}.
+If you need any custom linear warping (independent of the WCS, see @ref{Linear
warping basics}), you need to call the respective operation manually.
+These are described in @ref{Linear warps to be called explicitly}.
+Please note that you may not use both linear and non-linear modes
simultaneously.
+For example, you cannot scale or rotate the image while removing its
non-linear distortions at the same time.
+
+The following options are shared between both modes:
+
+@table @option
+@item --hstartwcs=INT
+Specify the first header keyword number (line) that should be used to read the
WCS information, see the full explanation in @ref{Invoking astcrop}.
+
+@item --hendwcs=INT
+Specify the last header keyword number (line) that should be used to read the
WCS information, see the full explanation in @ref{Invoking astcrop}.
+
+@item -C FLT
+@itemx --coveredfrac=FLT
+Depending on the warp, the output pixels that cover pixels on the edge of the
input image, or blank pixels in the input image, are not going to be fully
covered by input data.
+With this option, you can specify the acceptable covered fraction of such
pixels (any value between 0 and 1).
+If you only want output pixels that are fully covered by the input image area
(and are not blank), then you can set @option{--coveredfrac=1} (which is the
default!).
+Alternatively, a value of @code{0} will keep output pixels that are even
infinitesimally covered by the input.
+As a result, with @option{--coveredfrac=0}, the sum of the pixels in the input
and output images will be exactly the same.
+@end table
+
+@menu
+* Align pixels with WCS and remove distortions:: Default operation.
+* Linear warps to be called explicitly:: Other waprs.
+@end menu
+
+@node Align pixels with WCS and remove distortions, Linear warps to be called
explicitly, Invoking astwarp, Invoking astwarp
+@subsubsection Align pixels with WCS and remove distortions
+
+@cindex Resampling
+@cindex WCS distortion
+@cindex TPV distortion
+@cindex SIP distortion
+@cindex Non-linear distortion
+@cindex Align pixel and WCS coordinates
+When none of the linear warps@footnote{For linear warps, see @ref{Linear warps
to be called explicitly}.} are requested, Warp will align the input's pixel
axes with it's WCS axes.
+In the process, any possibly existing distortion is also removed (such as
@code{TPV} and @code{SIP}).
+Usually, the WCS axes are the Right Ascension and Declination in equatorial
coordinates.
+The output image's pixel grid is highly customizable through the options in
this section.
+
+@table @option
+@item -c FLT,FLT
+@itemx --center=FLT,FLT
+@cindex CRVALi
+@cindex Aligning an image
+WCS coordinates of the center of the central pixel of the output image.
+If @option{--center} isn't given, the output will have have the same central
WCS coordinate as the input.
+
+Usually, the WCS coordinates are Right Ascension and Declination (when the
first three characters of @code{CTYPE1} and @code{CTYPE2} are respectively
@code{RA-} and @code{DEC}).
+For more on the @code{CTYPEi} keyword values, see @code{--ctype}.
+
+@item -w INT,INT
+@itemx --widthinpix=INT,INT
+Pixel width and height of the output image.
+When @option{--widthinpix} isn't given, Warp will calculate the necessary size
of the output pixel grid to fully contain the input image.
+
+This option should be given two @emph{odd} integers that are greater than 1,
so that the output image can have a @emph{central} pixel.
+Recall that through the @option{--center} option, you specify the WCS
coordinate of the center of the central pixel.
+The central coordinate of an image with an even number of pixels will be on
the edge of two pixels, so a ``central'' pixel is not well defined.
+If any of the given values are even, Warp will automatically add a single
pixel (to make it an odd integer) and print a warning message.
-Warps must be specified as command-line options, either as (possibly multiple)
modular warpings (for example @option{--rotate}, or @option{--scale}), or
directly as a single raw matrix (with @option{--matrix}).
+@item -x FLT[,FLT]
+@itemx --cdelt=FLT[,FLT]
+@cindex CDELTi
+@cindex Pixel scale
+Coordinate deltas or increments (@code{CDELTi} in the FITS standard), or the
pixel scale in both dimensions.
+If a single value is given, it will be used for both axes.
+In this way, the output's pixels will be squares on the sky at the reference
point (as is usually expected!).
+When this option isn't given, Warp will read the input's pixel scale and
choose the larger of @code{CDELT1} or @code{CDELT2} so the output pixels are
square.
+
+Usually (when dealing with RA and Dec, and the @code{CUNITi}s have a value of
@code{deg}), the units of the given values are degrees/pixel.
+Warp allows you to easily convert from @emph{arcsec} to @emph{degrees} by
simply appending a @code{/3600} to the value.
+For example, for an output image of pixel scale @code{0.27} arcsec/pixel, you
can use @code{--cdelt=0.27/3600}.
+
+
+@item --ctype=STR,STR
+@cindex Align
+@cindex CTYPEi
+@cindex Resampling
+The coordinate types of the output (@code{CTYPE1} and @code{CTYPE2} keywords
in the FITS standard), separated by a comma.
+By default the value to this option is `@code{RA---TAN,DEC--TAN}'.
+Therefore if you don't call this option, the output WCS coordinates will be
Right Ascension and Declination, while the output's pojection will be
@url{https://en.wikipedia.org/wiki/Gnomonic_projection,Gnomonic}, also known as
Tantential (TAN).
+This combination is the most common in extragalactic imaging surveys.
+For other coordinates and projections in your output use other values, as
described below.
+
+According to the FITS standard version 4.0@footnote{FITS standard version 4.0:
@url{https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf}}:
@code{CTYPEi} is the
+``type for the Intermediate-coordinate Axis @mymath{i}.
+Any coordinate type that is not covered by this Standard or an officially
recognized FITS convention shall be taken to be linear.
+All non-linear coordinate system names must be expressed in `4–3' form: the
first four characters specify the coordinate type, the fifth character is a
hyphen (@code{-}), and the remaining three characters specify an algorithm code
for computing the world coordinate value.
+Coordinate types with names of fewer than four characters are padded on the
right with hyphens, and algorithm codes with fewer than three characters are
padded on the right with SPACE.
+Algorithm codes should be three characters''
+(see list of algorithm codes below).
+
+@cindex WCS Projections
+@cindex Projections (world coordinate system)
+You can use any of the projection algorithms (last three characters of each
coordinate's type) provided by your host WCSLIB (a mandatory dependency of
Gnuastro; see @ref{WCSLIB}).
+For a very elaborate and complete description of projection algorithms in the
FITS WCS standard, see @url{https://doi.org/10.1051/0004-6361:20021327,
Calabretta and Greisen 2002}.
+Wikipedia also has a nice article on
@url{https://en.wikipedia.org/wiki/Map_projection, Map projections}.
+As an example, WCSLIB 7.12 (released in September 2022) has the following
projection algorithms:
+
+@table @code
+@item AZP
+@cindex Zenithal/azimuthal projection
+Zenithal/azimuthal perspective
+@item SZP
+@cindex Slant zenithal projection
+Slant zenithal perspective
+@item TAN
+@cindex Gnomonic (tangential) projection
+Gnomonic (tangential)
+@item STG
+@cindex Stereographic projection
+Stereographic
+@item SIN
+@cindex Orthographic/synthesis projection
+Orthographic/synthesis
+@item ARC
+@cindex Zenithal/azimuthal equidistant projection
+Zenithal/azimuthal equidistant
+@item ZPN
+@cindex Zenithal/azimuthal polynomial projection
+Zenithal/azimuthal polynomial
+@item ZEA
+@cindex Zenithal/azimuthal equal area projection
+Zenithal/azimuthal equal area
+@item AIR
+@cindex Airy projection
+Airy
+@item CYP
+@cindex Cylindrical perspective projection
+Cylindrical perspective
+@item CEA
+@cindex Cylindrical equal area projection
+Cylindrical equal area
+@item CAR
+@cindex Plate carree projection
+Plate carree
+@item MER
+@cindex Mercator projection
+Mercator
+@item SFL
+@cindex Sanson-Flamsteed projection
+Sanson-Flamsteed
+@item PAR
+@cindex Parabolic projection
+Parabolic
+@item MOL
+@cindex Mollweide projection
+Mollweide
+@item AIT
+@cindex Hammer-Aitoff projection
+Hammer-Aitoff
+@item COP
+@cindex Conic perspective projection
+Conic perspective
+@item COE
+@cindex Conic equal area projection
+Conic equal area
+@item COD
+@cindex Conic equidistant projection
+Conic equidistant
+@item COO
+@cindex Conic orthomorphic projection
+Conic orthomorphic
+@item BON
+@cindex Bonne projection
+Bonne
+@item PCO
+@cindex Polyconic projection
+Polyconic
+@item TSC
+@cindex Tangential spherical cube projection
+Tangential spherical cube
+@item CSC
+@cindex COBE spherical cube projection
+COBE spherical cube
+@item QSC
+@cindex Quadrilateralized spherical cube projection
+Quadrilateralized spherical cube
+@item HPX
+@cindex HEALPix projection
+HEALPix
+@item XPH
+@cindex Butterfly projection
+@cindex HEALPix polar projection
+HEALPix polar, aka "butterfly"
+@end table
+
+
+
+@item --edgesampling=INT
+Number of extra samplings along the edge of a pixel.
+By default the value is @code{0} (the output pixel's polygon over the input
will be a quadrilateral (a polygon with four edges/vertices).
+
+Warp uses pixel mixing to derive the output pixel values, for a complete
introduction.
+See @ref{Resampling}, and in particular its later part on distortions.
+To account for this possible curvature due to distortion, you can use this
option.
+For example, @option{--edgesampling=1} will add one extra vertice in the
middle of each edge of the output pixel, producing an 8-vertice polygon.
+Similarly, @option{--edgesampling=5} will put 5 extra vertices along each
edge, thus sampling the shape (and possible curvature) of the output pixel over
an input pixel with @mymath{4+5\times4=24} vertice polygon.
+Since the polygon clipping will happen for every output pixel, a higher value
to this option can significantly reduce the running speed and increase the RAM
usage of Warp; so use it with caution: in most cases the default
@option{--edgesampling=0} is sufficient.
+@end table
+
+
+
+
+
+@node Linear warps to be called explicitly, , Align pixels with WCS and
remove distortions, Invoking astwarp
+@subsubsection Linear warps to be called explicitly
+
+Linear warps include operations like rotation, scaling, sheer and etc.
+For an introduction, see @ref{Linear warping basics}.
+These are warps that don't depend on the WCS of the image and should be
expliclitly requested.
+To align the input pixel coordinates with the WCS coordinates, see @ref{Align
pixels with WCS and remove distortions}.
+
+While they will correct any existing WCS based on the warp, they can also
operate on images without any WCS.
+For example, you have a mock image that doesn't (yet!) have its mock WCS, and
it has been created on an over-sampled grid and convolved with an over-sampled
PSF.
+In this scenario, you can use the @option{--scale} option to under-sample it
to your desired resolution.
+This is similar to the @ref{Sufi simulates a detection} tutorial.
+
+Linear warps must be specified as command-line options, either as (possibly
multiple) modular warpings (for example @option{--rotate}, or
@option{--scale}), or directly as a single raw matrix (with @option{--matrix}).
If specified together, the latter (direct matrix) will take precedence and all
the modular warpings will be ignored.
Any number of modular warpings can be specified on the command-line and
configuration files.
If more than one modular warping is given, all will be merged to create one
warping matrix.
@@ -19501,30 +19759,20 @@ So the coordinate center [0.0, 0.0] is half a pixel
away (in each axis) from the
The resampling that is done in Warp (see @ref{Resampling}) is done on the
coordinate axes and thus directly depends on the coordinate center.
In some situations this if fine, for example when rotating/aligning a real
image, all the edge pixels will be similarly affected.
But in other situations (for example when scaling an over-sampled mock image
to its intended resolution, this is not desired: you want the center of the
coordinates to be on the corner of the pixel.
-In such cases, you can use the @option{--centeroncorner} option which will
shift the center by @mymath{0.5} before the main warp, then shift it back by
@mymath{-0.5} after the main warp, see below.
-
+In such cases, you can use the @option{--centeroncorner} option which will
shift the center by @mymath{0.5} before the main warp, then shift it back by
@mymath{-0.5} after the main warp.
@table @option
-@item -a
-@itemx --align
-Align the image and celestial (WCS) axes given in the input.
-After it, the vertical image direction (when viewed in SAO DS9) corresponds
to the declination and the horizontal axis is the inverse of the Right
Ascension (RA).
-The inverse of the RA is chosen so the image can correspond to what you would
actually see on the sky and is common in most survey images.
-
-Align is internally treated just like a rotation (@option{--rotation}), but
uses the input image's WCS to find the rotation angle.
-Thus, if you have rotated the image before calling @option{--align}, you might
get unexpected results (because the rotation is defined on the original WCS).
-
@item -r FLT
@itemx --rotate=FLT
-Rotate the input image by the given angle in degrees: @mymath{\theta} in
@ref{Warping basics}.
+Rotate the input image by the given angle in degrees: @mymath{\theta} in
@ref{Linear warping basics}.
Note that commonly, the WCS structure of the image is set such that the RA is
the inverse of the image horizontal axis which increases towards the right in
the FITS standard and as viewed by SAO DS9.
So the default center for rotation is on the right of the image.
If you want to rotate about other points, you have to translate the warping
center first (with @option{--translate}) then apply your rotation and then
return the center back to the original position (with another call to
@option{--translate}, see @ref{Merging multiple warpings}.
@item -s FLT[,FLT]
@itemx --scale=FLT[,FLT]
-Scale the input image by the given factor(s): @mymath{M} and @mymath{N} in
@ref{Warping basics}.
+Scale the input image by the given factor(s): @mymath{M} and @mymath{N} in
@ref{Linear warping basics}.
If only one value is given, then both image axes will be scaled with the given
value.
When two values are given (separated by a comma), the first will be used to
scale the first axis and the second will be used for the second axis.
If you only need to scale one axis, use @option{1} for the axis you do not
need to scale.
@@ -19540,7 +19788,7 @@ Hence, if you want to flip by the second axis only, use
@option{--flip=0,1}.
@item -e FLT[,FLT]
@itemx --shear=FLT[,FLT]
-Shear the input image by the given value(s): @mymath{A} and @mymath{B} in
@ref{Warping basics}.
+Shear the input image by the given value(s): @mymath{A} and @mymath{B} in
@ref{Linear warping basics}.
If only one value is given, then both image axes will be sheared with the
given value.
When two values are given (separated by a comma), the first will be used to
shear the first axis and the second will be used for the second axis.
If you only need to shear along one axis, use @option{0} for the axis that
must be untouched.
@@ -19548,7 +19796,7 @@ The value(s) can also be written (on the command-line
or in configuration files)
@item -t FLT[,FLT]
@itemx --translate=FLT[,FLT]
-Translate (move the center of coordinates) the input image by the given
value(s): @mymath{c} and @mymath{f} in @ref{Warping basics}.
+Translate (move the center of coordinates) the input image by the given
value(s): @mymath{c} and @mymath{f} in @ref{Linear warping basics}.
If only one value is given, then both image axes will be translated by the
given value.
When two values are given (separated by a comma), the first will be used to
translate the first axis and the second will be used for the second axis.
If you only need to translate along one axis, use @option{0} for the axis that
must be untouched.
@@ -19556,7 +19804,7 @@ The value(s) can also be written (on the command-line
or in configuration files)
@item -p FLT[,FLT]
@itemx --project=FLT[,FLT]
-Apply a projection to the input image by the given values(s): @mymath{g} and
@mymath{h} in @ref{Warping basics}.
+Apply a projection to the input image by the given values(s): @mymath{g} and
@mymath{h} in @ref{Linear warping basics}.
If only one value is given, then projection will apply to both axes with the
given value.
When two values are given (separated by a comma), the first will be used to
project the first axis and the second will be used for the second axis.
If you only need to project along one axis, use @option{0} for the axis that
must be untouched.
@@ -19567,27 +19815,20 @@ The value(s) can also be written (on the command-line
or in configuration files)
The warp/transformation matrix.
All the elements in this matrix must be separated by commas(@key{,})
characters and as described above, you can also use fractions (a forward-slash
between two numbers).
The transformation matrix can be either a 2 by 2 (4 numbers), or a 3 by 3 (9
numbers) array.
-In the former case (if a 2 by 2 matrix is given), then it is put into a 3 by 3
matrix (see @ref{Warping basics}).
+In the former case (if a 2 by 2 matrix is given), then it is put into a 3 by 3
matrix (see @ref{Linear warping basics}).
@cindex NaN
The determinant of the matrix has to be non-zero and it must not contain any
non-number values (for example infinities or NaNs).
The elements of the matrix have to be written row by row.
-So for the general Homography matrix of @ref{Warping basics}, it should be
called with @command{--matrix=a,b,c,d,e,f,g,h,1}.
+So for the general Homography matrix of @ref{Linear warping basics}, it should
be called with @command{--matrix=a,b,c,d,e,f,g,h,1}.
The raw matrix takes precedence over all the modular warping options listed
above, so if it is called with any number of modular warps, the latter are
ignored.
-@item -c
-@itemx --centeroncorner
+@item --centeroncorner
Put the center of coordinates on the corner of the first (bottom-left when
viewed in SAO DS9) pixel.
This option is applied after the final warping matrix has been finalized:
either through modular warpings or the raw matrix.
See the explanation above for coordinates in the FITS standard to better
understand this option and when it should be used.
-@item --hstartwcs=INT
-Specify the first header keyword number (line) that should be used to read the
WCS information, see the full explanation in @ref{Invoking astcrop}.
-
-@item --hendwcs=INT
-Specify the last header keyword number (line) that should be used to read the
WCS information, see the full explanation in @ref{Invoking astcrop}.
-
@item -k
@itemx --keepwcs
@cindex WCSLIB
@@ -19595,14 +19836,6 @@ Specify the last header keyword number (line) that
should be used to read the WC
Do not correct the WCS information of the input image and save it untouched to
the output image.
By default the WCS (World Coordinate System) information of the input image is
going to be corrected in the output image so the objects in the image are at
the same WCS coordinates.
But in some cases it might be useful to keep it unchanged (for example to
correct alignments).
-
-@item -C FLT
-@itemx --coveredfrac=FLT
-Depending on the warp, the output pixels that cover pixels on the edge of the
input image, or blank pixels in the input image, are not going to be fully
covered by input data.
-With this option, you can specify the acceptable covered fraction of such
pixels (any value between 0 and 1).
-If you only want output pixels that are fully covered by the input image area
(and are not blank), then you can set @option{--coveredfrac=1}.
-Alternatively, a value of @code{0} will keep output pixels that are even
infinitesimally covered by the input(so the sum of the pixels in the input and
output images will be the same).
-
@end table
@@ -19655,7 +19888,7 @@ The Statistics program is designed for such situations.
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
-* Least squares fitting::
+* Least squares fitting:: Fitting with various parametric functions.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@end menu
@@ -24807,7 +25040,7 @@ Our aim is to put a radial profile of any functional
form @mymath{f(r)} over an
Hence we need to associate a radius/distance to every point in space.
Let's define the radial distance @mymath{r_{el}} as the distance on the major
axis to the center of an ellipse which is located at @mymath{i_c} and
@mymath{j_c} (in other words @mymath{r_{el}\equiv{a}}).
We want to find @mymath{r_{el}} of a point located at @mymath{(i,j)} (in the
image coordinate system) from the center of the ellipse with axis ratio
@mymath{q} and position angle @mymath{\theta}.
-First the coordinate system is rotated@footnote{Do not confuse the signs of
@mymath{sin} with the rotation matrix defined in @ref{Warping basics}.
+First the coordinate system is rotated@footnote{Do not confuse the signs of
@mymath{sin} with the rotation matrix defined in @ref{Linear warping basics}.
In that equation, the point is rotated, here the coordinates are rotated and
the point is fixed.} by @mymath{\theta} to get the new rotated coordinates of
that point @mymath{(i_r,j_r)}:
@dispmath{i_r(i,j)=+(i_c-i)\cos\theta+(j_c-j)\sin\theta}
@@ -24831,7 +25064,7 @@ In short, when a point is rotated in this order, we
first rotate it around the Z
Following the discussion in @ref{Merging multiple warpings}, we can define the
full rotation with the following matrix multiplication.
However, here we are rotating the coordinates, not the point.
Therefore, both the rotation angles and rotation order are reversed.
-We are also not using homogeneous coordinates (see @ref{Warping basics}) since
we are not concerned with translation in this context:
+We are also not using homogeneous coordinates (see @ref{Linear warping
basics}) since we are not concerned with translation in this context:
@dispmath{\left[\matrix{i_r\cr j_r\cr k_r}\right] =
\left[\matrix{cos\gamma&sin\gamma&0\cr -sin\gamma&cos\gamma&0\cr
0&0&1}\right]
@@ -29293,6 +29526,7 @@ If you use the Info version of this manual (see
@ref{Info}), you do not have to
* Labeled datasets:: Working with Segmented/labeled datasets.
* Convolution functions:: Library functions to do convolution.
* Interpolation:: Interpolate (over blank values possibly).
+* Warp library:: Warp pixel grid to a new one.
* Color functions:: Definitions and operations related to colors.
* Git wrappers:: Wrappers for functions in libgit2.
* Python interface:: Functions to help in writing Python wrappers.
@@ -33192,89 +33426,29 @@ The FITS standard defines the world coordinate system
(WCS) as a mechanism to as
For example, it can be used to convert pixel coordinates in an image to
celestial coordinates like the right ascension and declination.
The functions in this section are mainly just wrappers over CFITSIO, WCSLIB
and GSL library functions to help in common applications.
-In the two sub-sub-sections below, we will first explain the best approach to
avoid memory leaks when freeing a WCS structure.
-Afterwards, we'll describe all the WCS-related functions of Gnuastro.
-
-@menu
-* Freeing a WCS structure:: Demo on creating and properly freeing.
-* WCS Function and global constants:: All functions in gnuastro/wcs.h.
-@end menu
-
-@node Freeing a WCS structure, WCS Function and global constants, World
Coordinate System, World Coordinate System
-@subsubsection Freeing a WCS structure
-Through the functions in @ref{WCS Function and global constants} you can read
a WCS structure from a file (@code{gal_wcs_read}), an already-open CFITSIO FITS
pointer (@code{gal_wcs_read_fitsptr}) or create one yourself
(@code{gal_wcs_create}).
-In any of these cases, after you have completed your work with the WCS
structure, you should free it.
-
-To free a WCS pointer that comes from Gnuastro's functions, first free the WCS
structure's contents using WCSLIB's @code{wcsfree} function.
-Afterwards, free the actual structure with the @code{free} function of the
standard C library (@file{stdlib.h}).
-Remember that @code{wcsfree} only frees the contents, not the structure that
was allocated inside Gnuastro's library.
-
-Below you can see a fully working example program that will create a WCS
structure from scratch and will later free it properly (without any memory
leak).
-If you save the code below in a file called @file{wcs-create-free.c}, you can
compile and run it using @ref{BuildProgram} with this command:
@command{astbuildprog wcs-create-free.c}.
-
-@example
-@verbatim
-#include <stdio.h>
-#include <stdlib.h>
-#include <gnuastro/wcs.h>
-
-int
-main(void)
-{
- int status;
- size_t ndim=2;
- struct wcsprm *wcs;
- double crpix[]={50, 50};
- double pc[]={-1, 0, 0, 1};
- double cdelt[]={0.4, 0.4};
- double crval[]={178.23, 36.98};
- char *cunit[]={"deg", "deg"};
- char *ctype[]={"RA---TAN", "DEC--TAN"};
- int linearmatrix = GAL_WCS_LINEAR_MATRIX_PC;
-
- /* Allocate and fill the 'wcsprm' structure. */
- wcs=gal_wcs_create(crpix, crval, cdelt, pc, cunit,
- ctype, ndim, linearmatrix);
- printf("WCS structure created.\n");
-
- /* Free the contents (and check if there was a problem). */
- status=wcsfree(wcs);
- if(status)
- {
- /* Inform the user of the problem. */
- fprintf(stderr, "%s: WCSLIB wcsfree ERROR %d: %s",
- __func__, status, wcs_errmsg[status]);
-
- /* Return with a failure (don't continue with the
- rest of the function). */
- return EXIT_FAILURE;
- }
-
- /* Let the user know */
- printf("WCS structure's contents freed.\n");
-
- /* Free the actual 'wcsprm' structure. */
- free(wcs);
- printf("WCS structure freed.\n");
-
- /* Return successfully. */
- return EXIT_SUCCESS;
-}
-@end verbatim
-@end example
-
-@node WCS Function and global constants, , Freeing a WCS structure, World
Coordinate System
-@subsubsection WCS Function and macros
+@cindex Thread safety
+@cindex WCSLIB thread safety
+[@strong{Tread safety}] Since WCSLIB version 5.18 (released in January 2018),
most WCSLIB functions are thread
safe@footnote{@url{https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/threads.html}}.
+Gnuastro has high-level functions to easily spin-off threads and speed up your
programs.
+For a fully working example see @ref{Library demo - multi-threaded operation}.
+However you still need to be cautious in the following scenarios below.
+@itemize
+@item
+Many users or operating systems may still use an older version.
+@item
+The @code{wcsprm} structure of WCSLIB is not thread-safe: you can't use the
same pointer on multiple threads.
+For example if you use @code{gal_wcs_img_to_world} simultaneously on multiple
threads, you shouldn't pass the same @code{wcsprm} structure pointer.
+You can use @code{gal_wcs_copy} to keep and use separate copies the main
structure within each thread, and later free the copies with
@code{gal_wcs_free}.
+@end itemize
The full set of functions and global constants that are defined by Gnuastro's
@file{gnuastro/wcs.h} are described below.
-For freeing a WCS structure after reading or creating one, see @ref{Freeing a
WCS structure}.
-
-@deffn Macro GAL_WCS_DISTORTION_TPD
-@deffnx Macro GAL_WCS_DISTORTION_SIP
-@deffnx Macro GAL_WCS_DISTORTION_TPV
-@deffnx Macro GAL_WCS_DISTORTION_DSS
-@deffnx Macro GAL_WCS_DISTORTION_WAT
-@deffnx Macro GAL_WCS_DISTORTION_INVALID
+
+@deffn {Global integer} GAL_WCS_DISTORTION_TPD
+@deffnx {Global integer} GAL_WCS_DISTORTION_SIP
+@deffnx {Global integer} GAL_WCS_DISTORTION_TPV
+@deffnx {Global integer} GAL_WCS_DISTORTION_DSS
+@deffnx {Global integer} GAL_WCS_DISTORTION_WAT
+@deffnx {Global integer} GAL_WCS_DISTORTION_INVALID
@cindex WCS distortion
@cindex Distortion, WCS
@cindex TPD WCS distortion
@@ -33290,13 +33464,13 @@ TPD is a superset of all these, hence it has both
prior and sequeal distortion c
More information is given in the documentation of @code{dis.h}, from the
WCSLIB
manual@footnote{@url{https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/dis_8h.html}}.
@end deffn
-@deffn Macro GAL_WCS_COORDSYS_EQB1950
-@deffnx Macro GAL_WCS_COORDSYS_EQJ2000
-@deffnx Macro GAL_WCS_COORDSYS_ECB1950
-@deffnx Macro GAL_WCS_COORDSYS_ECJ2000
-@deffnx Macro GAL_WCS_COORDSYS_GALACTIC
-@deffnx Macro GAL_WCS_COORDSYS_SUPERGALACTIC
-@deffnx Macro GAL_WCS_COORDSYS_INVALID
+@deffn {Global integer} GAL_WCS_COORDSYS_EQB1950
+@deffnx {Global integer} GAL_WCS_COORDSYS_EQJ2000
+@deffnx {Global integer} GAL_WCS_COORDSYS_ECB1950
+@deffnx {Global integer} GAL_WCS_COORDSYS_ECJ2000
+@deffnx {Global integer} GAL_WCS_COORDSYS_GALACTIC
+@deffnx {Global integer} GAL_WCS_COORDSYS_SUPERGALACTIC
+@deffnx {Global integer} GAL_WCS_COORDSYS_INVALID
@cindex Galactic coordinate system
@cindex Ecliptic coordinate system
@cindex Equatorial coordinate system
@@ -33310,9 +33484,9 @@ Recognized WCS coordinate systems in Gnuastro.
In the equatorial and ecliptic coordinates, @code{B1950} stands for the
Besselian 1950 epoch and @code{J2000} stands for the Julian 2000 epoch.
@end deffn
-@deffn Macro GAL_WCS_LINEAR_MATRIX_PC
-@deffnx Macro GAL_WCS_LINEAR_MATRIX_CD
-@deffnx Macro GAL_WCS_LINEAR_MATRIX_INVALID
+@deffn {Global integer} GAL_WCS_LINEAR_MATRIX_PC
+@deffnx {Global integer} GAL_WCS_LINEAR_MATRIX_CD
+@deffnx {Global integer} GAL_WCS_LINEAR_MATRIX_INVALID
Identifiers of the linear transformation matrix: either in the @code{PCi_j} or
the @code{CDi_j} formalism.
For more, see the description of @option{--wcslinearmatrix} in @ref{Input
output options}.
@end deffn
@@ -33322,11 +33496,59 @@ For more, see the description of
@option{--wcslinearmatrix} in @ref{Input output
Limit of rounding for floating point errors.
@end deffn
+@deftypefun {struct wcsprm *} gal_wcs_create (double @code{*crpix}, double
@code{*crval}, double @code{*cdelt}, double @code{*pc}, char @code{**cunit},
char @code{**ctype}, size_t @code{ndim}, int @code{linearmatrix})
+Given all the most common standard components of the WCS standard, construct a
@code{struct wcsprm}, initialize and set it for future processing.
+See the FITS WCS standard for more on these keywords.
+All the arrays must have @code{ndim} elements with them except for @code{pc}
which should have @code{ndim*ndim} elements (a square matrix).
+Also, @code{cunit} and @code{ctype} are arrays of strings.
+If @code{GAL_WCS_LINEAR_MATRIX_CD} is passed to @code{linearmatrix} then the
output WCS structure will have a CD matrix (even though you have given a PC and
CDELT matrix as input to this function).
+Otherwise, the output will have a PC and CDELT matrix (which is the
recommended format by WCSLIB).
+
+@example
+@verbatim
+#include <stdio.h>
+#include <stdlib.h>
+#include <gnuastro/wcs.h>
+
+int
+main(void)
+{
+ int status;
+ size_t ndim=2;
+ struct wcsprm *wcs;
+ double crpix[]={50, 50};
+ double pc[]={-1, 0, 0, 1};
+ double cdelt[]={0.4, 0.4};
+ double crval[]={178.23, 36.98};
+ char *cunit[]={"deg", "deg"};
+ char *ctype[]={"RA---TAN", "DEC--TAN"};
+ int linearmatrix = GAL_WCS_LINEAR_MATRIX_PC;
+
+ /* Allocate and fill the 'wcsprm' structure. */
+ wcs=gal_wcs_create(crpix, crval, cdelt, pc, cunit,
+ ctype, ndim, linearmatrix);
+ printf("WCS structure created.\n");
+
+ /*... Add any operation with the WCS structure here ...*/
+
+ /* Free the WCS structure. */
+ gal_wcs_free(wcs);
+ printf("WCS structure freed.\n");
+
+ /* Return successfully. */
+ return EXIT_SUCCESS;
+}
+@end verbatim
+@end example
+@end deftypefun
+
@deftypefun {struct wcsprm *} gal_wcs_read_fitsptr (fitsfile @code{*fptr}, int
@code{linearmatrix}, size_t @code{hstartwcs}, size_t @code{hendwcs}, int
@code{*nwcs})
-[@strong{Not thread-safe}] Return the WCSLIB @code{wcsprm} structure that is
read from the CFITSIO @code{fptr} pointer to an opened FITS file.
+Return the WCSLIB @code{wcsprm} structure that is read from the CFITSIO
@code{fptr} pointer to an opened FITS file.
+With older WCSLIB versions (in particular below version 5.18) this function
may not be thread-safe.
+
Also put the number of coordinate representations found into the space that
@code{nwcs} points to.
To read the WCS structure directly from a filename, see @code{gal_wcs_read}
below.
-After processing has finished, you should free the WCS structure, for doing
that properly, see @ref{Freeing a WCS structure}.
+After processing has finished, you should free the WCS structure that this
function returns with @code{gal_wcs_free}.
The @code{linearmatrix} argument takes one of three values: @code{0},
@code{GAL_WCS_LINEAR_MATRIX_PC} and @code{GAL_WCS_LINEAR_MATRIX_CD}.
It will determine the format of the WCS when it is later written to file with
@code{gal_wcs_write} or @code{gal_wcs_write_in_fitsptr} (which is called by
@code{gal_fits_img_write})
@@ -33349,18 +33571,15 @@ Therefore, be sure to not call this function
simultaneously (over multiple threa
Also put the number of coordinate representations found into the space that
@code{nwcs} points to.
Please see @code{gal_wcs_read_fitsptr} for more.
-After processing has finished, you should free the WCS structure, for doing
that properly, see @ref{Freeing a WCS structure}.
+After processing has finished, you should free the WCS structure that this
function returns with @code{gal_wcs_free}.
@end deftypefun
-@deftypefun {struct wcsprm *} gal_wcs_create (double @code{*crpix}, double
@code{*crval}, double @code{*cdelt}, double @code{*pc}, char @code{**cunit},
char @code{**ctype}, size_t @code{ndim}, int @code{linearmatrix})
-Given all the most common standard components of the WCS standard, construct a
@code{struct wcsprm}, initialize and set it for future processing.
-See the FITS WCS standard for more on these keywords.
-All the arrays must have @code{ndim} elements with them except for @code{pc}
which should have @code{ndim*ndim} elements (a square matrix).
-Also, @code{cunit} and @code{ctype} are arrays of strings.
-If @code{GAL_WCS_LINEAR_MATRIX_CD} is passed to @code{linearmatrix} then the
output WCS structure will have a CD matrix (even though you have given a PC and
CDELT matrix as input to this function).
-Otherwise, the output will have a PC and CDELT matrix (which is the
recommended format by WCSLIB).
-
-For a fully working example of creating a WCS structure with this function and
how to later free it, see @ref{Freeing a WCS structure}.
+@deftypefun void gal_wcs_free (struct wcsprm @code{*wcs})
+Free the contents @emph{and} the space that @code{wcs} points to.
+WCSLIB's @code{wcsfree} function only frees the contents of the @code{wcsprm}
structure, not the actual pointer.
+However, Gnuastro's @code{wcsprm} creation and reading functions allocate the
structure also.
+This higher-level function therefore simplifies the processing.
+A complete working example is given in the description of
@code{gal_wcs_create}.
@end deftypefun
@deftypefun {char *} gal_wcs_dimension_name (struct wcsprm @code{*wcs}, size_t
@code{dimension})
@@ -34678,7 +34897,7 @@ The input (@code{in}) is an array containing the
coordinates (two values) of eac
So @code{in} should have @code{2*n} elements.
The output (@code{ordinds}) is an array with @code{n} elements specifying the
indices in order.
This array must have been allocated before calling this function.
-The indexes are output for more generic usage, for example in a homographic
transform (necessary in warping an image, see @ref{Warping basics}), the
necessary order of vertices is the same for all the pixels.
+The indexes are output for more generic usage, for example in a homographic
transform (necessary in warping an image, see @ref{Linear warping basics}), the
necessary order of vertices is the same for all the pixels.
In other words, only the positions of the vertices change, not the way they
need to be ordered.
Therefore, this function would only be necessary once.
@@ -34830,7 +35049,7 @@ functions, a global variable or structure are also
necessary as described
below.
@deffn {Global variable} {gal_qsort_index_single}
-@cindex Thread-safety
+@cindex Thread safety
@cindex Multi-threaded operation
Pointer to an array (for example @code{float *} or @code{int *}) to use as
a reference in @code{gal_qsort_index_single_TYPE_d} or
@@ -36396,7 +36615,7 @@ channel borders. Other pixels in the image will not be
touched. Hence, it
is much faster.
@end deftypefun
-@node Interpolation, Color functions, Convolution functions, Gnuastro library
+@node Interpolation, Warp library, Convolution functions, Gnuastro library
@subsection Interpolation (@file{interpolate.h})
@cindex Sky line
@@ -36637,7 +36856,146 @@ blank. To see if any blank (non-interpolated)
elements remain, you can use
@code{gal_blank_present} on @code{in} after this function is finished.
@end deftypefun
-@node Color functions, Git wrappers, Interpolation, Gnuastro library
+
+
+
+
+@cindex Warp
+@cindex Align
+@cindex Resampling
+@cindex WCS distortion
+@cindex Non-linear distortion
+@node Warp library, Color functions, Interpolation, Gnuastro library
+@subsection Warp library (@file{warp.h})
+
+Warping an image to a new pixel grid is commonly necessary as part of
astronomical data redution, for an introduction, see @ref{Warp}.
+For details of how we resample the old pixel grid to the new pixel grid, see
@ref{Resampling}.
+Gnuastro's Warp program uses the following functions for its default mode
(when no linear warps are requested).
+Through the following functions, you can directly access those features in
your own custom programs.
+The linear warping operations of the Warp program aren't yet brought into the
library.
+If you need them please get in touch with us at @code{bug-gnuastro@@gnu.org}.
+For a demo please refer to @ref{Library demo - resampling an image to a new
grid}.
+
+You are free to provide any valid WCS keywords to the functions defined in
this library using the @code{gal_warp_wcsalign_t} data type.
+This might be used to align the input image to the standard WCS grid,
potentially changing the pixel scale, removing any valid WCS non-linear
distortion available, and projecting to any valid WCS projection type.
+Further details of the warp library functions and parameters are shown below:
+
+@deftp {Type (C @code{struct})} gal_warp_wcsalign_t
+The main data container for inputs, output and internal variables to simplify
the WCS-aligning functions.
+Due to the large number of input variables, this structure makes it easy to
call the main functions.
+Similar to @code{gal_data_t}, the @code{gal_warp_wcsalign_t} is a structure
@code{typedef}'d as a new type, see @ref{Library data container}.
+Please note that this structure has elements that are @emph{allocated}
dynamically and must be freed after usage.
+@code{gal_warp_wcsalign_free} only frees the internal variables, so you are
responsible for freeing your own inputs (like @code{cdelt}, @code{input} and
etc) and the output.
+The internal variables are cached here to cut cpu-intensive computations.
+The structure and each of its elements are defined below:
+
+@example
+typedef struct
+@{
+ /* Arguments given (and later freed) by the caller. */
+ gal_data_t *ctype;
+ gal_data_t *cdelt;
+ gal_data_t *input;
+ gal_data_t *widthinpix;
+ gal_data_t *center;
+ size_t numthreads;
+ double coveredfrac;
+ size_t edgesampling;
+
+ /* Output (must be freed by caller) */
+ gal_data_t *output;
+
+ /* Internal variables (allocated and freed internally) */
+ size_t v0;
+ size_t nhor;
+ size_t ncrn;
+ size_t gcrn;
+ int isccw;
+ gal_data_t *vertices;
+@} gal_warp_wcsalign_t;
+@end example
+
+@table @code
+@item gal_data_t *ctype
+The output's projection type.
+The dataset has to have the type @code{GAL_TYPE_STRING}, containing exactly
two strings.
+Both strings will be directly passed to WCSLIB and should conform to the FITS
standard's @code{CTYPEi} keywords, see the description of @option{--ctype} in
@ref{Align pixels with WCS and remove distortions}.
+For example, @code{"RA---TAN"} and @code{"DEC--TAN"}, or @code{"RA---HPX"} and
@code{"DEC--HPX"}.
+
+@item gal_data_t *cdelt
+Output pixel scale (size of pixel in the WCS units: value to @code{CUNITi}
keywords in FITS, usually degrees).
+The dataset should have a type of @code{GAL_TYPE_FLOAT64} and contain exactly
two values.
+Hint: to convert arcsec to degrees, just divide by 3600.
+
+@item gal_data_t *input
+The input dataset.
+This dataset must contain both the image array of type
@code{GAL_TYPE_FLOAT64}, and @code{input->wcs} should not be @code{NULL} for
the WCS-aligning operations to work, see @ref{Library demo - resampling an
image to a new grid}.
+
+@item gal_data_t *widthinpix
+Output image size (width and height) in number of pixels.
+If a @code{NULL} pointer is passed, the WCS-aligning operations will estimate
the output image size internally such that it contains the full input.
+This dataset should have a type of @code{GAL_TYPE_SIZE_T} and contain exactly
two @emph{odd} values.
+This ensures that the center of the central pixel lies at the requested
central coordinate (note that an image with an even number of pixels doesn't
have a ``central'' pixel!
+
+@item gal_data_t *center
+WCS coordinate of the center of the central pixel of the output.
+The units depend on the WCS, for example if the @code{CUNITi} keywords are
@code{deg}, it is in degrees.
+This dataset should have a type of @code{GAL_TYPE_FLOAT64} and contain exactly
two values.
+
+@item size_t numthreads
+Number of threads to use during the WCS aligning operations.
+If the given value is @code{0}, the library will calculate the number of
available threads at run-time.
+The @code{warp} library functions are @emph{thread-safe} so you can freely
enjoy the merits of parallel processing.
+
+@item double coveredfrac
+Acceptable fraction of output pixel that is covered by input pixels.
+The value should be between 0 and 1 (inclusive).
+If the area of an output pixel is covered by less than this fraction, its
value will be @code{NaN}.
+For more, see the description of @option{--coveredfrac} in @ref{Invoking
astwarp}.
+
+@item size_t edgesampling
+Set the number of extra vertices along each edge of the output pixel's polygon
to account for potential curvature due to projection or distortion.
+A value of @code{0} is usually enough for this (so the pixel is only defined
by a four vertice polygon.
+Greater values increase memory usage and program execution time.
+For more, plese see the description of @option{--edgesampling} in @ref{Align
pixels with WCS and remove distortions}.
+@end table
+@end deftp
+
+
+@deftypefun void gal_warp_wcsalign (gal_warp_wcsalign_t *wa)
+A high-level function to align the input dataset's pixels to its WCS
coordinates and write the result in @code{wa->output}.
+This function assumes that the input variables have already been set in the
@code{wa} structure.
+The input variables are clearly shown in the definition of
@code{gal_warp_wcsalign_t}.
+It will call the lower level functions below to do the job and will free the
internal variables afterwards.
+@end deftypefun
+
+The following low-level functions are called from the high-level
@code{gal_warp_wcsalign} function.
+They are provided here in scenarios where fine grain control over the thread
workflow is necessary, see @ref{Multithreaded programming}.
+
+@deftypefun void gal_warp_wcsalign_init (gal_warp_wcsalign_t *wa)
+Low-level function to initialize all the elements inside the @code{wa}
structure assuming that the input variables have been set.
+The input variables are clearly shown in the definition of
@code{gal_warp_wcsalign_t}.
+This includes sanity checking the input arguments, as well as allocating the
output image's empty pixels (that can be filled with
@code{gal_warp_wcsalign_onpix}, possibly on threads).
+@end deftypefun
+
+@deftypefun void gal_warp_wcsalign_onpix (gal_warp_wcsalign_t *nl, size_t ind)
+Low-level function that fills pixel @code{ind} (counting from 0) in the
already initialized output image.
+@end deftypefun
+
+@deftypefun {void *} gal_warp_wcsalign_onthread (void *inparam)
+Low-level worker function that can be passed to the high-level
@code{gal_threads_spin_off} or the lower-level @code{pthread_create} with some
modifications, see @ref{Multithreaded programming}.
+@end deftypefun
+
+@deftypefun void gal_warp_wcsalign_free (gal_warp_wcsalign_t *wa)
+Low-level function to free the internal variables inside @code{wa} only.
+The caller must free the input pointers themselves, this function will not
free them (they may be necessary in other parts of the caller's higher-level
architecture).
+@end deftypefun
+
+
+
+
+
+@node Color functions, Git wrappers, Warp library, Gnuastro library
@subsection Color functions (@file{color.h})
@cindex Colors
@@ -37159,6 +37517,7 @@ use those that are generated by Gnuastro after
@command{make check} in the
* Library demo - inspecting neighbors:: Inspect the neighbors of a pixel.
* Library demo - multi-threaded operation:: Doing an operation on threads.
* Library demo - reading and writing table columns:: Simple Column I/O.
+* Library demo - resampling an image to a new grid::
@end menu
@node Library demo - reading a image, Library demo - inspecting neighbors,
Library demo programs, Library demo programs
@@ -37406,7 +37765,7 @@ main(void)
-@node Library demo - reading and writing table columns, , Library demo -
multi-threaded operation, Library demo programs
+@node Library demo - reading and writing table columns, Library demo -
resampling an image to a new grid, Library demo - multi-threaded operation,
Library demo programs
@subsection Library demo - reading and writing table columns
Tables are some of the most common inputs to, and outputs of programs. This
@@ -37446,6 +37805,7 @@ the column before writing your program (see
@ref{Numeric data types}), so
type checking and copying to a specific type will not be necessary.
@example
+@verbatim
#include <stdio.h>
#include <stdlib.h>
@@ -37453,7 +37813,7 @@ type checking and copying to a specific type will not
be necessary.
int
main(void)
-@{
+{
/* File names and column names (which may also be numbers). */
char *c1_name="NAME1", *c2_name="NAME2", *c3_name="10";
char *inname="input.fits", *hdu="1", *outname="out.fits";
@@ -37484,7 +37844,7 @@ main(void)
counter=1;
for(col=columns; col!=NULL; col=col->next)
switch(counter++)
- @{
+ {
case 1: /* First column: we want it as int32_t. */
c1=gal_data_copy_to_new_type(col, GAL_TYPE_INT32);
array1 = c1->array;
@@ -37497,25 +37857,25 @@ main(void)
case 3: /* Third column: it MUST be double. */
if(col->type!=GAL_TYPE_FLOAT64)
- @{
+ {
fprintf(stderr, "Column %s must be float64 type, it is "
"%s", c3_name, gal_type_name(col->type, 1));
exit(EXIT_FAILURE);
- @}
+ }
array3 = col->array;
break;
- @}
+ }
/* As an example application we will just print them out. In the
* meantime (just for a simple demonstration), change the first
* array value to the counter and multiply the second by 10. */
for(i=0;i<c1->size;++i)
- @{
+ {
printf("%zu: %d + %f + %f = %f\n", i+1, array1[i], array2[i],
array3[i], array1[i]+array2[i]+array3[i]);
array1[i] = i+1;
array2[i] *= 10;
- @}
+ }
/* Link the first two columns as a list. */
c1->next = c2;
@@ -37537,12 +37897,110 @@ main(void)
gal_list_data_free(columns);
gal_list_str_free(column_ids, 0); /* strings were not allocated. */
return EXIT_SUCCESS;
-@}
+}
+@end verbatim
@end example
+@node Library demo - resampling an image to a new grid, , Library demo -
reading and writing table columns, Library demo programs
+@subsection Library demo - resampling an image to a new grid
+
+Gnuastro's warp library allows you to resample an image from a grid to another
entirely using the WCSLIB.
+You may use it to remove non-linear distortions from an image, change pixel
scale, projection type, or just align the image to get it ready for
@emph{stacking} several images.
+Here you can find an example program that demonstrates the usage of the
respective functions in @ref{Warp library}.
+
+To compile the program, put the code below in a plain-text file (let's assume
its called @file{align.c}) and use @ref{BuildProgram} with this command:
`@command{astbuildprog align.c}'.
+
+@example
+@verbatim
+#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";
+
+ /* Output file name. */
+ char *outname="output.fits";
+
+ /* RA/Dec of the center of the central pixel of output. Please
+ change the center based on your input. */
+ double center[]={52.55401335, -28.74252337};
+
+ /* Coordinate and Projection algorithms of output. */
+ char *ctype[2]={"RA---TAN", "DEC--TAN"};
+
+ /* Output pixel scale (in units of degrees/pixel). */
+ double cdelt[]={0.27/3600, 0.27/3600};
+
+ /* For intermediate steps. */
+ size_t two=2;
+ gal_warp_wcsalign_t nl;
+
+ /* Set the width (and height!) of the output in pixels (as a 1D and
+ 2 element 'gal_data_t'). When it is NULL, the library will
+ calculate the appropriate width to fully fit the input image
+ after alignment. */
+ nl.widthinpix=NULL;
+
+ /* Set the number of threads to use. If the value is '0', the
+ library will estimate the maximum available threads at
+ run-time on the host operating system. */
+ nl.numthreads=0;
+
+
+ /* Read the input image and its WCS. */
+ nl.input=gal_array_read_one_ch_to_type(filename, hdu, NULL,
+ GAL_TYPE_FLOAT64, -1, 0);
+ nl.input->wcs=gal_wcs_read(filename, hdu, 0, 0, 0, &nl.input->nwcs);
+
+
+ /* Prepare the warp input structure. */
+ nl.coveredfrac=1; nl.edgesampling=0;
+ nl.ctype=gal_data_alloc(ctype, GAL_TYPE_STRING, 1, &two, NULL, 1,
+ -1, 0, NULL, NULL, NULL);
+ nl.cdelt=gal_data_alloc(cdelt, GAL_TYPE_FLOAT64, 1, &two, NULL, 1,
+ -1, 0, NULL, NULL, NULL);
+ nl.center=gal_data_alloc(center, GAL_TYPE_FLOAT64, 1, &two, NULL, 1,
+ -1, 0, NULL, NULL, NULL);
+
+
+ /* Do the warp, then convert it to a 32-bit float. */
+ gal_warp_wcsalign(&nl);
+ nl.output=gal_data_copy_to_new_type_free(nl.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(nl.output, outname, NULL, "warp-demo");
+
+
+ /* Remove the pointers to arrays that we didn't allocate (and thus,
+ should not be freed by 'gal_data_free' below). */
+ nl.cdelt->array=nl.center->array=nl.ctype->array=NULL;
+
+
+ /* Clean up and return. */
+ gal_data_free(nl.cdelt); gal_data_free(nl.ctype);
+ gal_data_free(nl.input); gal_data_free(nl.output);
+ gal_data_free(nl.center); gal_data_free(nl.widthinpix);
+
+ /* Give control back to the operating system. */
+ return EXIT_SUCCESS;
+}
+@end verbatim
+@end example
+
+
+
@@ -39435,26 +39893,24 @@ This will enable you to start your branch (work) from
the most recent commit and
@node Gnuastro programs list, Other useful software, Developing, Top
@appendix Gnuastro programs list
-GNU Astronomy Utilities @value{VERSION}, contains the following
-programs. They are sorted in alphabetical order and a short description is
-provided for each program. The description starts with the executable names
-in @file{thisfont} followed by a pointer to the respective section in
-parenthesis. Throughout this book, they are ordered based on their context,
-please see the top-level contents for contextual ordering (based on what
-they do).
+GNU Astronomy Utilities @value{VERSION}, contains the following programs.
+They are sorted in alphabetical order and a short description is provided for
each program.
+The description starts with the executable names in @file{thisfont} followed
by a pointer to the respective section in parenthesis.
+Throughout this book, they are ordered based on their context, please see the
top-level contents for contextual ordering (based on what they do).
@table @asis
@item Arithmetic
(@file{astarithmetic}, see @ref{Arithmetic}) For arithmetic operations on
multiple (theoretically unlimited) number of datasets (images).
-It has a large and growing set of arithmetic, mathematical, and even
statistical operators (for example @command{+}, @command{-}, @command{*},
@command{/}, @command{sqrt}, @command{log}, @command{min}, @command{average},
@command{median}).
+It has a large and growing set of arithmetic, mathematical, and even
statistical operators (for example @command{+}, @command{-}, @command{*},
@command{/}, @command{sqrt}, @command{log}, @command{min}, @command{average},
@command{median}, see @ref{Arithmetic operators}).
@item BuildProgram
-(@file{astbuildprog}, see @ref{BuildProgram}) Compile, link and run programs
that depend on the Gnuastro library (see @ref{Gnuastro library}).
+(@file{astbuildprog}, see @ref{BuildProgram}) Compile, link and run custom C
programs that depend on the Gnuastro library (see @ref{Gnuastro library}).
This program will automatically link with the libraries that Gnuastro depends
on, so there is no need to explicitly mention them every time you are compiling
a Gnuastro library dependent program.
@item ConvertType
(@file{astconvertt}, see @ref{ConvertType}) Convert astronomical data files
(FITS or IMH) to and from several other standard image and data formats, for
example TXT, JPEG, EPS or PDF.
+Optionally, it is also possible to add vector graphics markers over the output
image (for example circles from catalogs containing RA or Dec).
@item Convolve
(@file{astconvolve}, see @ref{Convolve}) Convolve (blur or smooth) data with a
given kernel in spatial and frequency domain on multiple threads.
@@ -39464,8 +39920,8 @@ Convolve can also do de-convolution to find the
appropriate kernel to PSF-match
(@file{astcosmiccal}, see @ref{CosmicCalculator}) Do cosmological
calculations, for example the luminosity distance, distance modulus, comoving
volume and many more.
@item Crop
-(@file{astcrop}, see @ref{Crop}) Crop region(s) from an image and stitch
several images if necessary.
-Inputs can be in pixel coordinates or world coordinates.
+(@file{astcrop}, see @ref{Crop}) Crop region(s) from one or many image(s) and
stitch several images if necessary.
+Input coordinates can be in pixel coordinates or world coordinates.
@item Fits
(@file{astfits}, see @ref{Fits}) View and manipulate FITS file extensions and
header keywords.
@@ -39497,14 +39953,16 @@ It uses a technique to detect very faint and diffuse,
irregularly shaped signal
@item Statistics
(@file{aststatistics}, see @ref{Statistics}) Statistical calculations on the
input dataset (column in a table, image or datacube).
+This includes man operations like generating histogram, sigma clipping, least
squares fitting and etc.
@item Table
-(@file{asttable}, @ref{Table}) Convert FITS binary and ASCII tables into other
such tables, print them on the command-line, save them in a plain text file, or
get the FITS table information.
+(@file{asttable}, @ref{Table}) Convert FITS binary and ASCII tables into other
such tables, print them on the command-line, save them in a plain text file, do
arithmetic on the columns or get the FITS table information.
+For a full list of operations, see @ref{Operation precedence in Table}.
@item Warp
(@file{astwarp}, see @ref{Warp}) Warp image to new pixel grid.
-Any projective transformation or Homography can be applied to the input images.
-
+By default it will align the pixel and WCS coordinates, removing any
non-linear WCS distortions.
+Any linear warp (projective transformation or Homography) can also be applied
to the input images by explicitly calling the respective operation.
@end table
The programs listed above are designed to be highly modular and generic.
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 440b2f88..168771a5 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -144,12 +144,14 @@ libgnuastro_la_SOURCES = \
txt.c \
type.c \
units.c \
+ warp.c \
wcs.c
+
# Installed headers, note that we are not blindly including all '.h' files
# in the $(headersdir) directory. Some of the header files don't need to be
# installed.
@@ -190,6 +192,7 @@ pkginclude_HEADERS = gnuastro/config.h \
$(headersdir)/txt.h \
$(headersdir)/type.h \
$(headersdir)/units.h \
+ $(headersdir)/warp.h \
$(headersdir)/wcs.h
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index b575b039..275b7b27 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
must be included before the C++ preparations below */
#include <gnuastro/data.h>
#include <gnuastro/blank.h>
+#include <gnuastro/polygon.h>
/* C++ Preparations */
#undef __BEGIN_C_DECLS
@@ -72,6 +73,28 @@ gal_dimension_num_neighbors(size_t ndim);
#define GAL_DIMENSION_FLT_TO_INT(FLT) ( (FLT)-(long)(FLT) > 0.5f \
? (long)(FLT)+1 : (long)(FLT) )
+/* A pixel's center is an integer value. This function will give the
+ integer value that is nearest to a floating point number. Works for both
+ positive and negative values and includes floating point errors.
+
+ WARP_NEARESTINT_HALFHIGHER(0.5f) --> 1.0f
+*/
+#define GAL_DIMENSION_NEARESTINT_HALFHIGHER(D) \
+ ( ceil( (D) ) - (D) > 0.5f + GAL_POLYGON_ROUND_ERR ? \
+ ceil((D))-1.0f : ceil((D)) )
+
+/* Similar to 'GAL_DIMENSION_NEARESTINT_HALFHIGHER' but:
+
+ GAL_DIMENSION_NEARESTINT_HALFLOWER(0.5f) --> 0.0f; */
+#define GAL_DIMENSION_NEARESTINT_HALFLOWER(D) \
+ ( ceil( (D) ) - (D) > 0.5f - GAL_POLYGON_ROUND_ERR ? \
+ ceil((D))-1.0f : ceil((D)) )
+
+#define GAL_DIMENSION_CEILWITHERR(D) \
+ ( fabs( nearbyint((D)) - (D) ) < GAL_POLYGON_ROUND_ERR ? \
+ nearbyint((D)) : nearbyint((D))+1 )
+
+
void
gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
diff --git a/lib/gnuastro/warp.h b/lib/gnuastro/warp.h
new file mode 100644
index 00000000..e736b1f5
--- /dev/null
+++ b/lib/gnuastro/warp.h
@@ -0,0 +1,109 @@
+/*********************************************************************
+Warp -- Warp pixels of one dataset to another pixel grid.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Corresponding author:
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
+Contributing author(s):
+ Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Copyright (C) 2022 Free Software Foundation, Inc.
+
+Gnuastro 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.
+
+Gnuastro 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 Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#ifndef __GAL_WARP_H__
+#define __GAL_WARP_H__
+
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+#include <wcslib/wcs.h>
+
+#include <gnuastro/data.h>
+
+/* C++ Preparations */
+#undef __BEGIN_C_DECLS
+#undef __END_C_DECLS
+#ifdef __cplusplus
+# define __BEGIN_C_DECLS extern "C" {
+# define __END_C_DECLS }
+#else
+# define __BEGIN_C_DECLS /* empty */
+# define __END_C_DECLS /* empty */
+#endif
+/* End of C++ preparations */
+
+
+
+
+
+/* Actual header contants (the above were for the Pre-processor). */
+__BEGIN_C_DECLS /* From C++ preparations */
+
+
+typedef struct
+{
+ /* Arguments given (and later freed) by the caller. */
+ gal_data_t *ctype; /* Type of the coordinates. */
+ gal_data_t *cdelt; /* Pixel scale of the output image. */
+ gal_data_t *input; /* Pointer to input data structure. */
+ gal_data_t *center; /* Center of the image in RA and Dec. */
+ gal_data_t *widthinpix; /* Output image width and height in pixels. */
+ size_t numthreads; /* Number of threads to use. */
+ double coveredfrac; /* Acceptable fraction of output covered. */
+ size_t edgesampling; /* Order of samplings along each pixel edge. */
+
+ /* Output (must be freed by caller) */
+ gal_data_t *output; /* Pointer to output data structure. */
+
+ /* Internal variables (allocated and freed internally) */
+ size_t v0; /* The first vertical corner in each pixel. */
+ size_t nhor; /* No. of vertices on top side of output. */
+ size_t ncrn; /* No. of total vertices around each pixel. */
+ size_t gcrn; /* Gap between corners of each row. */
+ int isccw; /* Rotation orientation of pixel edges. */
+ gal_data_t *vertices; /* Stores all vertice coords of output img. */
+} gal_warp_wcsalign_t;
+
+
+
+
+
+/* Create the empty output WCS-ready image */
+void
+gal_warp_wcsalign_init(gal_warp_wcsalign_t *wa);
+
+
+/* Fill nonlinear output by pixel */
+void
+gal_warp_wcsalign_onpix(gal_warp_wcsalign_t *wa, size_t ind);
+
+
+void *
+gal_warp_wcsalign_onthread(void *inparam);
+
+
+/* Spin-off the threads and finalize the output 'gal_data_t' image in
+ 'nl->output' */
+void
+gal_warp_wcsalign(gal_warp_wcsalign_t *wa);
+
+
+/* Clean up the entire struct */
+void
+gal_warp_wcsalign_free(gal_warp_wcsalign_t *wa);
+
+
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif /* __GAL_WARP_H__ */
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 27c69c7b..6523f83c 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -112,6 +112,9 @@ gal_wcs_create(double *crpix, double *crval, double *cdelt,
double *pc, char **cunit, char **ctype, size_t ndim,
int linearmatrix);
+void
+gal_wcs_free(struct wcsprm *wcs);
+
char *
gal_wcs_dimension_name(struct wcsprm *wcs, size_t dimension);
diff --git a/lib/warp.c b/lib/warp.c
new file mode 100644
index 00000000..b77714d9
--- /dev/null
+++ b/lib/warp.c
@@ -0,0 +1,970 @@
+/*********************************************************************
+Warp -- Warp pixels of one dataset to another pixel grid.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Corresponding author:
+ Pedram Ashofteh Ardakani <pedramardakani@pm.me>
+Contributing author(s):
+ Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Copyright (C) 2022 Free Software Foundation, Inc.
+
+Gnuastro 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.
+
+Gnuastro 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 Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <errno.h>
+#include <error.h>
+#include <stdio.h>
+
+#include <gnuastro/wcs.h>
+#include <gnuastro/type.h>
+#include <gnuastro/warp.h>
+#include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
+#include <gnuastro/polygon.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
+
+
+
+
+
+/* Macros */
+#define WARP_NEXT_ODD(D) \
+ ( (size_t)( ceil((D)) ) + ( (size_t)( ceil((D)) )%2==0 ? 1 : 0 ) )
+
+
+#define WARP_WCSALIGN_H(IND,ES,IS1) \
+ (size_t)( ( (IND)%(IS1) ) * ( (ES)+1 ) \
+ + ( (IND)/(IS1) ) * (1+(IS1)*( (ES)+1 )) )
+
+
+#define WARP_WCSALIGN_V0(ES,IS0,IS1) \
+ (size_t)( 1+(IS0)+(IS1)*( (IS0)+1 )*( (ES)+1) )
+
+
+#define WARP_WCSALIGN_V(IND,ES,V0,IS1) \
+ (size_t)( (V0)+(ES)*( (IND)+(IND)/(IS1) ) )
+
+
+
+
+
+/* Generate the points on the outer boundary of a dsize[0] x dsize[1]
+ matrix and return the array */
+static gal_data_t *
+warp_alloc_perimeter(gal_data_t *input)
+{
+ /* Low level variables */
+ size_t ind, i;
+ gal_data_t *pcrn;
+ double *x=NULL, *y=NULL;
+ size_t is0=input->dsize[0];
+ size_t is1=input->dsize[1];
+ int quietmmap=input->quietmmap;
+ size_t minmapsize=input->minmapsize;
+
+ /* High level variables */
+ size_t npcrn=2*(is0+is1);
+
+ pcrn=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &npcrn, NULL, 0,
+ minmapsize, quietmmap, NULL, NULL, NULL);
+ pcrn->next=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &npcrn, NULL, 0,
+ minmapsize, quietmmap, NULL, NULL, NULL);
+
+ /* Find outermost pixel coordinates of the input image. Cover two
+ corners at once to shorten the loop */
+ x=pcrn->array; y=pcrn->next->array;
+
+ /* Top and bottom */
+ ind=0;
+ for(i=is1+1; i--;)
+ {
+ x[ind]=i+0.5f; y[ind]=0.5f; ind++;
+ x[ind]=i+0.5f; y[ind]=is0+0.5f; ind++;
+ }
+
+ /* Left and right */
+ for(i=is0-1; i--;)
+ {
+ x[ind]=0.5f; y[ind]=1.5f+i; ind++;
+ x[ind]=0.5f+is1; y[ind]=1.5f+i; ind++;
+ }
+
+ /* Sanity check: let's make sure we have correctly covered the input
+ perimeter */
+ if(ind!=npcrn)
+ error(EXIT_FAILURE, 0, "%s: the input img perimeter of size <%zu> "
+ "is not covered correctly. Currently on ind <%zu>.", __func__,
+ npcrn, ind);
+
+ return pcrn;
+}
+
+
+
+
+
+/* Create a base image with WCS consisting of the basic geometry
+ keywords */
+static void
+warp_wcsalign_init_output(gal_warp_wcsalign_t *wa)
+{
+ /* Low level variables */
+ size_t i, nkcoords, *osize;
+ gal_data_t *input=wa->input;
+ int quietmmap=input->quietmmap;
+ size_t minmapsize=input->minmapsize;
+ struct wcsprm *bwcs=NULL, *rwcs=NULL;
+ gal_data_t *kcoords=NULL, *pcrn=NULL, *converted=NULL;
+ double ocrpix[2], *xkcoords, *ykcoords, *x=NULL, *y=NULL;
+ double pmin[2]={DBL_MAX, DBL_MAX}, pmax[2]={-DBL_MAX, -DBL_MAX};
+
+ /* Base WCS default parameters */
+ double pc[4]={-1, 0, 0, 1};
+ double rcrpix[2]={1.0, 1.0};
+ char **ctype=wa->ctype->array;
+ struct wcsprm *iwcs=input->wcs;
+ double *cdelt=wa->cdelt->array;
+ double *center=wa->center->array;
+
+ /* High level variables */
+ char *cunit[2]={iwcs->cunit[0], iwcs->cunit[1]};
+
+ /* Determine the output image size */
+ size_t iminr=GAL_BLANK_SIZE_T, imaxr=GAL_BLANK_SIZE_T; /* indexs of */
+ size_t imind=GAL_BLANK_SIZE_T, imaxd=GAL_BLANK_SIZE_T; /* extreme-um */
+ double pminr=DBL_MAX, pmind=DBL_MAX, pmaxr=-DBL_MAX, pmaxd=-DBL_MAX;
+
+ /* Create the reference WCS */
+ rwcs=gal_wcs_create(rcrpix, center, cdelt, pc, cunit, ctype, 2,
+ GAL_WCS_LINEAR_MATRIX_PC);
+
+ /* Calculate the outer boundary of the input. */
+ pcrn=warp_alloc_perimeter(input);
+ converted=gal_wcs_img_to_world(pcrn, iwcs, 0);
+
+ /* Get the minimum/maximum of the outer boundary. */
+ x=converted->array; y=converted->next->array;
+ for(i=converted->size; i--;)
+ {
+ if(x[i]<pminr) { pminr=x[i]; iminr=i; }
+ if(y[i]<pmind) { pmind=y[i]; imind=i; }
+ if(x[i]>pmaxr) { pmaxr=x[i]; imaxr=i; }
+ if(y[i]>pmaxd) { pmaxd=y[i]; imaxd=i; }
+ }
+
+ /* Prepare the key world coorinates and change to image coordinates
+ later. We are doing this to determine the CRPIX and NAXISi size for
+ the final image */
+ nkcoords=5;
+ kcoords=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nkcoords, NULL, 0,
+ minmapsize, quietmmap, NULL, NULL, NULL);
+ kcoords->next=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nkcoords, NULL,
+ 0, minmapsize, quietmmap, NULL, NULL, NULL);
+ xkcoords=kcoords->array; ykcoords=kcoords->next->array;
+ xkcoords[0]=x[iminr]; ykcoords[0]=y[iminr]; /* min RA */
+ xkcoords[1]=x[imaxr]; ykcoords[1]=y[imaxr]; /* max RA */
+ xkcoords[2]=x[imind]; ykcoords[2]=y[imind]; /* min Dec */
+ xkcoords[3]=x[imaxd]; ykcoords[3]=y[imaxd]; /* max Dec */
+ xkcoords[4]=center[0]; ykcoords[4]=center[1]; /* Image center */
+
+ /* Convert to pixel coords */
+ gal_wcs_world_to_img(kcoords, rwcs, 1);
+
+ /* Determine output image size */
+ if( wa->widthinpix )
+ osize=wa->widthinpix->array;
+ else
+ {
+ /* Automatic: the first four coordinates are the extreme-um RA/Dec */
+ for(i=4; i--;)
+ {
+ pmin[0] = xkcoords[i] < pmin[0] ? xkcoords[i] : pmin[0];
+ pmin[1] = ykcoords[i] < pmin[1] ? ykcoords[i] : pmin[1];
+ pmax[0] = xkcoords[i] > pmax[0] ? xkcoords[i] : pmax[0];
+ pmax[1] = ykcoords[i] > pmax[1] ? ykcoords[i] : pmax[1];
+ }
+
+ /* Size must be odd so the image would have a center value. Also, the
+ indices are swapped since number of columns defines the horizontal
+ part of the center and vice versa. */
+ osize=gal_pointer_allocate(GAL_TYPE_SIZE_T, 2, 0,
+ __func__, "osize");
+ osize[0] = WARP_NEXT_ODD(pmax[1]-pmin[1]);
+ osize[1] = WARP_NEXT_ODD(pmax[0]-pmin[0]);
+ }
+
+ /* Set the CRPIX value
+
+ Note: os1 is number of columns, so we use it to define CRPIX in the
+ horizontal axis, and vice versa */
+ ocrpix[0]= 1.5f + osize[1]/2.0f - xkcoords[4];
+ ocrpix[1]= 1.5f + osize[0]/2.0f - ykcoords[4];
+
+ /* Create the base WCS */
+ bwcs=gal_wcs_create(ocrpix, center, cdelt, pc, cunit,
+ ctype, 2, GAL_WCS_LINEAR_MATRIX_PC);
+
+ /* Make sure that the size is reasonable (i.e., less than 100000 pixels
+ on a side). This can happen when a wrong central coordinate is
+ requested. */
+ if(osize[0]>100000 || osize[1]>100000)
+ error(EXIT_SUCCESS, 0, "%s: the output image size (%zu x %zu pixels) "
+ "is unreasonably large. This may be due to a mistake in the "
+ "given central coordinate compared to the input image (the "
+ "given center is too far from the image)", __func__, osize[1],
+ osize[0]);
+
+ /* Create the output image dataset with the base WCS */
+ wa->output=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 2, osize, bwcs, 0,
+ minmapsize, quietmmap, "Aligned", NULL, NULL);
+
+ /* Free */
+ wcsfree(bwcs);
+ wcsfree(rwcs);
+ gal_list_data_free(pcrn);
+ gal_list_data_free(kcoords);
+ gal_list_data_free(converted);
+ if(wa->widthinpix==NULL) free(osize);
+
+ /* Must be freed after wcsfree() is called! */
+ free(bwcs);
+ free(rwcs);
+}
+
+
+
+
+
+static gal_data_t *
+warp_wcsalign_init_vertices(gal_warp_wcsalign_t *wa)
+{
+ size_t es=wa->edgesampling;
+
+ size_t ind, i, j, ix, iy;
+ gal_data_t *vertices=NULL;
+ double gap=1.0f/(es+1.0f);
+ size_t os0=wa->output->dsize[0];
+ size_t os1=wa->output->dsize[1];
+ double *x=NULL, *y=NULL, row, col;
+ uint8_t quietmmap=wa->input->quietmmap;
+ size_t minmapsize=wa->input->minmapsize;
+
+ size_t nvcrn=es*os0;
+ size_t nhcrn=es*os1+os1+1;
+ size_t v0=WARP_WCSALIGN_V0(es,os0,os1);
+ size_t nvertices=nvcrn*(os1+1)+nhcrn*(os0+1);
+
+ /* Now create all subpixels based on the edgesampling '-E' option */
+ vertices=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nvertices, NULL, 0,
+ minmapsize, quietmmap, "OutputRA", NULL, NULL);
+ vertices->next=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nvertices,
+ NULL, 0, minmapsize, quietmmap,
+ "OutputDec", NULL, NULL);
+ x=vertices->array; y=vertices->next->array;
+
+ for(ind=os0*os1; ind--;)
+ {
+ row=ind%os1;
+ col=floor(ind/os1);
+ ix=WARP_WCSALIGN_H(ind,es,os1);
+ iy=WARP_WCSALIGN_V(ind,es,v0,os1);
+
+ /* Bottom left */
+ x[ix]= 0.5f+row;
+ y[ix]= 0.5f+col;
+
+ for(i=es; i--;)
+ {
+ /* Horizontal */
+ j=ix+i+1;
+ x[j]=0.5f+row+gap+i*gap;
+ y[j]=0.5f+col;
+
+ /* Vertical */
+ j=iy+i;
+ x[j]=0.5f+row;
+ y[j]=0.5f+col+gap+i*gap;
+ }
+ }
+
+ /* Top */
+ for(i=nhcrn; i--;)
+ {
+ j=v0-nhcrn+i;
+ x[j]=0.5f+gap*i;
+ y[j]=0.5f+os0;
+ }
+
+ /* Right */
+ for(ind=os1-1; ind<os1*os0; ind+=os1)
+ {
+ row=(size_t)(ind%os1); col=(size_t)(ind/os1);
+
+ iy=WARP_WCSALIGN_V(ind,es,v0,os1);
+ ix=WARP_WCSALIGN_H(ind,es,os1);
+
+ /* Bottom right */
+ j=ix+es+1;
+ x[j]=0.5f+os1;
+ y[j]=0.5f+col;
+
+ /* Right vertice */
+ for(i=es; i--;)
+ {
+ j=iy+es+i;
+ x[j]=0.5f+os1;
+ y[j]=0.5f+col+gap+i*gap;
+ }
+ }
+
+ return vertices;
+}
+
+
+
+
+
+static void
+warp_check_output_orientation(gal_warp_wcsalign_t *wa)
+{
+ size_t gcrn=wa->gcrn;
+ size_t es=wa->edgesampling;
+ double *vx=wa->vertices->array;
+ double *vy=wa->vertices->next->array;
+ size_t indices[4]={ 0, es+1, gcrn+es+1, gcrn };
+ double *temp=gal_pointer_allocate(GAL_TYPE_FLOAT64, 8, 0, __func__, NULL);
+
+ temp[ 0 ]=vx[ indices[0] ]; temp[ 1 ]=vy[ indices[0] ];
+ temp[ 2 ]=vx[ indices[1] ]; temp[ 3 ]=vy[ indices[1] ];
+ temp[ 4 ]=vx[ indices[2] ]; temp[ 5 ]=vy[ indices[2] ];
+ temp[ 6 ]=vx[ indices[3] ]; temp[ 7 ]=vy[ indices[3] ];
+
+ wa->isccw=gal_polygon_is_counterclockwise(temp, 4);
+
+ /* Clean up */
+ free(temp);
+}
+
+
+
+
+
+static double *
+warp_pixel_perimeter_ccw(gal_warp_wcsalign_t *wa, size_t ind)
+{
+ /* Low-level variables */
+ size_t i, j, hor, ver, ic;
+ double *xcrn=NULL, *ycrn=NULL, *ocrn=NULL;
+
+ /* High-level variables */
+ size_t v0=wa->v0;
+ size_t gcrn=wa->gcrn;
+ size_t ncrn=wa->ncrn;
+ size_t es=wa->edgesampling;
+ size_t os1=wa->output->dsize[1];
+
+ /* Set ocrn, the corners of each output pixel */
+ xcrn=wa->vertices->array;
+ ycrn=wa->vertices->next->array;
+ ocrn=gal_pointer_allocate(GAL_TYPE_FLOAT64, (2*ncrn), 0, __func__,
+ "ocrn");
+
+ /* Index of surrounding vertices for this pixel */
+ hor=WARP_WCSALIGN_H(ind, es, os1);
+ ver=WARP_WCSALIGN_V(ind, es, v0, os1);
+
+ /* All four corners
+
+ WARNING: this block of code highly depends on the ordering, take
+ extra care while refactoring
+
+ ocrn: all output edges transformed into the input image pixel
+ coordiantes
+
+ ic: index (position in array) of the current pixel edges
+
+ io: index of the output (reference) pixel edges
+
+
+ left edge -> +---------+ <- top edge
+ | |
+ | |
+ | |
+ bottom edge -> +---------+ <- right edge */
+
+
+ /* bottom left */
+ i=0;
+ j=hor;
+ ocrn[2*i]=xcrn[j]; ocrn[2*i+1]=ycrn[j];
+
+ /* bottom right */
+ i=es+1;
+ j=hor+es+1;
+ ocrn[2*i]=xcrn[j]; ocrn[2*i+1]=ycrn[j];
+
+ /* top right */
+ i=2*(es+1);
+ j=hor+es+1+gcrn;
+ ocrn[2*i]=xcrn[j]; ocrn[2*i+1]=ycrn[j];
+
+ /* top left */
+ i=3*(es+1);
+ j=hor+gcrn;
+ ocrn[2*i]=xcrn[j]; ocrn[2*i+1]=ycrn[j];
+
+ /* Sampling corners of the output pixel on the input image */
+ for(i=es; i--;)
+ {
+ /* bottom vertice: 0*(es+1)+(i+1) */
+ ic=i+1;
+ j=hor+i+1;
+ ocrn[2*ic]=xcrn[j]; ocrn[2*ic+1]=ycrn[j];
+
+ /* right vertice: 1*(es+1)+(i+1) */
+ ic=i+2+es;
+ j=ver+es+i;
+ ocrn[2*ic]=xcrn[j]; ocrn[2*ic+1]=ycrn[j];
+
+ /* top vertice: 2*(es+1)+(i+1) */
+ ic=i+3+2*es;
+ j=hor+es+gcrn-i;
+ ocrn[2*ic]=xcrn[j]; ocrn[2*ic+1]=ycrn[j];
+
+ /* left vertice: 3*(es+1)+(i+1) */
+ ic=i+4+3*es;
+ j=ver+es-i-1;
+ ocrn[2*ic]=xcrn[j]; ocrn[2*ic+1]=ycrn[j];
+ }
+
+ return ocrn;
+}
+
+
+
+
+
+static double *
+warp_pixel_perimeter_cw(gal_warp_wcsalign_t *wa, size_t ind)
+{
+ size_t i, hor, ver, ic;
+ double *xcrn=NULL, *ycrn=NULL, *ocrn=NULL;
+
+ size_t gcrn=wa->gcrn;
+ size_t ncrn=wa->ncrn;
+ size_t es=wa->edgesampling;
+ size_t os1=wa->output->dsize[1];
+
+ /* High level variables */
+ size_t v0=wa->v0;
+
+ /* Set ocrn, the corners of each output pixel */
+ xcrn=wa->vertices->array;
+ ycrn=wa->vertices->next->array;
+ ocrn=gal_pointer_allocate(GAL_TYPE_FLOAT64, 2*ncrn, 0, __func__, "ocrn");
+
+ /* Index of surrounding vertices for this pixel */
+ hor=WARP_WCSALIGN_H(ind, es, os1);
+ ver=WARP_WCSALIGN_V(ind, es, v0, os1);
+
+ /* All four corners: same as the counter clockwise method */
+
+ /* top left <- previously bottom left */
+ i=0;
+ ocrn[ 2*i ]=xcrn[ hor+gcrn ]; /* xcrn[ hor ] */
+ ocrn[ 2*i+1 ]=ycrn[ hor+gcrn ]; /* ycrn[ hor ] */
+
+ /* top right <- previously bottom right */
+ i=es+1;
+ ocrn[ 2*i ]=xcrn[ hor+es+1+gcrn ]; /* xcrn[ hor+es+1 ] */
+ ocrn[ 2*i+1 ]=ycrn[ hor+es+1+gcrn ]; /* ycrn[ hor+es+1 ] */
+
+ /* bottom right <- previously top right */
+ i=2*(es+1);
+ ocrn[ 2*i ]=xcrn[ hor+es+1 ]; /* xcrn[ hor+es+1+gcrn ] */
+ ocrn[ 2*i+1 ]=ycrn[ hor+es+1 ]; /* ycrn[ hor+es+1+gcrn ] */
+
+ /* bottom left <- previously top left */
+ i=3*(es+1);
+ ocrn[ 2*i ]=xcrn[ hor ]; /* xcrn=[ hor+gcrn ] */
+ ocrn[ 2*i+1 ]=ycrn[ hor ]; /* ycrn=[ hor+gcrn ] */
+
+ /* Sampling corners of the output pixel on the input image */
+ for(i=es; i--;)
+ {
+ /* top vertice 0*(es+1)+(i+1) <- previously bottom left */
+ ic=i+1;
+ ocrn[ 2*ic ]=xcrn[ hor+i+1+gcrn ]; /* xcrn=[ hor+i+1 ] */
+ ocrn[ 2*ic+1 ]=ycrn[ hor+i+1+gcrn ]; /* ycrn=[ hor+i+1 ] */
+
+ /* right vertice 1*(es+1)+(i+1) <- previously bottom right */
+ ic=i+2+es;
+ ocrn[ 2*ic ]=xcrn[ ver+2*es-1-i ]; /* xcrn[ ver+es+i ] */
+ ocrn[ 2*ic+1 ]=ycrn[ ver+2*es-1-i ]; /* ycrn[ ver+es+i ] */
+
+ /* bottom vertice 2*(es+1)+(i+1) <- previously top right */
+ ic=i+3+2*es;
+ ocrn[ 2*ic ]=xcrn[ hor+es-i ]; /* xcrn[ hor+es+gcrn-i ] */
+ ocrn[ 2*ic+1 ]=ycrn[ hor+es-i ]; /* ycrn[ hor+es+gcrn-i ] */
+
+ /* left vertice 3*(es+1)+(i+1) <- previously top left */
+ ic=i+4+3*es;
+ ocrn[ 2*ic ]=xcrn[ ver+i ]; /* xcrn[ ver+es-i-1 ] */
+ ocrn[ 2*ic+1 ]=ycrn[ ver+i ]; /* ycrn[ ver+es-i-1 ] */
+ }
+
+ return ocrn;
+}
+
+
+
+
+
+static void
+warp_wcsalign_check_2d(gal_data_t *in, uint8_t type, const char *func,
+ char *name, char *comment)
+{
+ if(!in)
+ error(EXIT_FAILURE, 0, "%s: no '%s' specified. %s", func,
+ name, comment);
+ if(in->size!=2)
+ error(EXIT_FAILURE, 0, "%s: '%s' takes exactly 2 values, currently "
+ "detected %zu values", func, name, in->size);
+ if(in->type!=type)
+ error(EXIT_FAILURE, 0, "%s: '%s' must have a type of '%s' but has "
+ "type '%s'", func, name, gal_type_name(type, 1),
+ gal_type_name(in->type, 1));
+}
+
+
+
+
+
+static void
+warp_wcsalign_init_params(gal_warp_wcsalign_t *wa, const char *func)
+{
+ size_t *tmp=NULL;
+
+ /* Check if input and 'wa' are not NULL! */
+ if(wa==NULL) error(EXIT_FAILURE, 0, "%s: 'wa' structure is NULL", func);
+ if(wa->input==NULL) error(EXIT_FAILURE, 0, "%s: input is NULL", func);
+
+ /* This function assumes the input is double precision. */
+ if(wa->input->type != GAL_TYPE_FLOAT64)
+ error(EXIT_FAILURE, 0, "%s: input must have a double precision "
+ "floating point type, but its type is '%s', you can use "
+ "'gal_data_copy_to_new_type' or "
+ "'gal_data_copy_to_new_type_free' for the conversion", func,
+ gal_type_name(wa->input->type, 1));
+
+ /* Check the 2D inputs. */
+ warp_wcsalign_check_2d(wa->ctype, GAL_TYPE_STRING, func, "ctype",
+ "Any pair of valid WCSLIB ctype is "
+ "allowed, e.g. 'RA---TAN, DEC--TAN'");
+ warp_wcsalign_check_2d(wa->cdelt, GAL_TYPE_FLOAT64, func, "cdelt",
+ "This is the pixel scale in degrees");
+ warp_wcsalign_check_2d(wa->center, GAL_TYPE_FLOAT64, func, "center",
+ "This is the output image center in degrees");
+
+ /* check 'widthinpix', it can be null for automatic detection */
+ if(wa->widthinpix)
+ {
+ warp_wcsalign_check_2d(wa->widthinpix, GAL_TYPE_SIZE_T, func,
+ "widthinpix", "This is the output "
+ "image size");
+ tmp=wa->widthinpix->array;
+ if(tmp[0]%2==0 || tmp[1]%2==0)
+ error(EXIT_FAILURE, 0, "%s: 'widthinpix' takes exactly 2 ODD "
+ "values, detected EVEN value in %zux%zu", func, tmp[0],
+ tmp[1]);
+ }
+
+ /* Check 'coveredfrac'. */
+ if(isnan( wa->coveredfrac ))
+ error(EXIT_FAILURE, 0, "%s: no 'coveredfrac' specified. This is the "
+ "acceptable fraction of output covered", func);
+ if(wa->coveredfrac<0.0f || wa->coveredfrac>1.0f)
+ error(EXIT_FAILURE, 0, "%s: coveredfrac takes exactly on positive "
+ "value less than or equal to 1.0, but it is given a value "
+ "of %f", func, wa->coveredfrac);
+
+ /* Check 'edgesampling', can't compare to '0' since it has meaning, can't
+ check if negative since it is an unsigned type */
+ if(wa->edgesampling==GAL_BLANK_SIZE_T)
+ error(EXIT_FAILURE, 0, "%s: no 'edgesampling' specified. This is the "
+ "Order of samplings along each pixel edge", func);
+ if(wa->edgesampling>999)
+ error(EXIT_FAILURE, 0, "%s: edgesampling takes zero OR a positive "
+ "integer value of type 'size_t', <%zu> is too big which might "
+ "be a bad cast", func, wa->edgesampling);
+
+ /* If 'numthreads' is 0, use the number of threads available to the
+ system. */
+ if(wa->numthreads==GAL_BLANK_SIZE_T || wa->numthreads==0)
+ wa->numthreads=gal_threads_number();
+
+ /* Initialize the internal parameters */
+ wa->vertices=NULL;
+ wa->isccw=GAL_BLANK_INT;
+ wa->v0=GAL_BLANK_SIZE_T;
+ wa->nhor=GAL_BLANK_SIZE_T;
+ wa->ncrn=GAL_BLANK_SIZE_T;
+ wa->gcrn=GAL_BLANK_SIZE_T;
+}
+
+
+
+
+
+/* Convert the necessary vertice coordinates. */
+static void *
+warp_wcsalign_init_convert(void *in_prm)
+{
+ /* Low-level definitions to be done first. */
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ gal_warp_wcsalign_t *wa = (gal_warp_wcsalign_t *)tprm->params;
+
+ /* Higher-level variables. */
+ gal_data_t *vertices=NULL;
+ double *xarr=wa->vertices->array;
+ int quietmmap=wa->vertices->quietmmap;
+ double *yarr=wa->vertices->next->array;
+ size_t minmapsize=wa->vertices->minmapsize;
+ size_t first, size, nt=wa->numthreads, vsize=wa->vertices->size;
+
+ /* WCSLIB's conversion functions write intermediate processing steps in
+ the 'wcsprm', so each thread should use its own copy. */
+ struct wcsprm *iwcs=gal_wcs_copy(wa->input->wcs);
+ struct wcsprm *owcs=gal_wcs_copy(wa->output->wcs);
+
+ /* Find the first vertice index to use in this thread. For the last
+ thread, the size will not be pre-defined. */
+ size = vsize/nt;
+ first = vsize/nt*tprm->id;
+ if(tprm->id==nt-1 && nt>1) size=vsize-(nt-1)*size;
+
+ /* For a check:
+ printf("%s: thread-%zu: %zu, %zu\n", __func__,
+ tprm->id, first, size);
+ */
+
+ /* Allocate the non-allocated vertices table for this thread. */
+ gal_list_data_add_alloc(&vertices, xarr+first, GAL_TYPE_FLOAT64,
+ 1, &size, NULL, 0, minmapsize, quietmmap,
+ NULL, NULL, NULL);
+ gal_list_data_add_alloc(&vertices, yarr+first, GAL_TYPE_FLOAT64,
+ 1, &size, NULL, 0, minmapsize, quietmmap,
+ NULL, NULL, NULL);
+ gal_list_data_reverse(&vertices); /* '_add' is last-in-first-out. */
+
+ /* Convert the coordinates. */
+ gal_wcs_img_to_world(vertices, owcs, 1);
+ gal_wcs_world_to_img(vertices, iwcs, 1);
+
+ /* Clean up: since the 'array' pointer is within a larger allocated
+ array, we shouldn't free it when freeing the table, so we'll set it to
+ NULL. */
+ vertices->array=vertices->next->array=NULL;
+ gal_list_data_free(vertices);
+ gal_wcs_free(iwcs);
+ gal_wcs_free(owcs);
+
+ /* Wait for all the other threads to finish, then return. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+/* Determine the final image size and allocate the output array
+ accordingly.
+
+ 'is0' indicates number of rows available in the input fits image,
+ while 'is1' indicates nummber of columns. The same goes for 'os0'
+ and 'os1'. See the following figure:
+
+ +------------------------+
+ / | |
+ | | N |
+ | | ^ |
+ | | | |
+ is0 < | | |
+ | | E <-----+ |
+ | | |
+ | | input image |
+ \ | |
+ +------------------------+
+ \__________ ___________/
+ v
+ is1
+
+ NOTE: please keep in mind that dsize[1] is NAXIS1 and dsize[0] is
+ NAXIS2 in FITS format.
+
+ In preparations, 'pcrn' is of type 'gal_data_t' linked lists. The
+ first list holds RA coords, and the next is Dec coords. This
+ variable is filled with the outer-most pixel coordinates. The
+ purpose of this variable is to hold the min and max RA and Dec
+ coordinates 'temporarily', so we can determine the output image size
+ later.
+
+ nhcrn and nvcrn:
+
+ 'nhcrn' stands for number of horizontal corners, and similarly
+ 'nvcrn' stands for number of vertical corners. Note that number of
+ 'corners' is different from number of 'pixels'. Each pixel has many
+ corners.
+
+ Also bear in mind, to keep from counting repeated corners on the
+ image edges, we let the horizontal corners devour the first and
+ last vertical corners. See the following figure:
+
+ hc6 hc7 hc8 hc9 hc10
+ +-----+-----+-----+-----+
+ / | |
+ | | img: 4x3 |
+ | x vc2 x vc4
+ is0=3 | | N |
+ | | ^ |
+ | x vc1 | x vc3
+ | | E <---+ |
+ \ | |
+ +-----+-----+-----+-----+
+ hc1 hc2 hc3 hc4 hc5
+
+ \_____________________/
+ is1=4
+
+ In the example above, for an image of 4x3, there will be 5
+ horizontal and 2 vertical corners in each axis, hence we would have
+ 10 horizontal and 4 vertical corners: 'nhcrn=2*(is1+1)' and
+ 'nvcrn=2*(is0-1)'. The total number of corners will be:
+ 'nhcrn+nvcrn=2*(is0-1)+2*(is1+1)=2*(is0+is1)=2*(4+3)=14'.
+
+
+ After finding out the min and max RA and Decs, the 'pcrn' is
+ projected back to pixel coordinates.
+*/
+void
+gal_warp_wcsalign_init(gal_warp_wcsalign_t *wa)
+{
+ size_t os0, os1, gcrn;
+ gal_data_t *input=wa ? wa->input : NULL;
+ size_t es = wa ? wa->edgesampling : GAL_BLANK_SIZE_T;
+
+ /* Run a sanity check on the input parameters */
+ warp_wcsalign_init_params(wa, __func__);
+
+ /* Create the output image */
+ warp_wcsalign_init_output(wa);
+ os0=wa->output->dsize[0];
+ os1=wa->output->dsize[1];
+ gcrn=1+os1*(es+1);
+
+ /* Set up the output image corners in pixel coords */
+ wa->vertices=warp_wcsalign_init_vertices(wa);
+
+ /* Project the output image corners to the input image pixel coords. We
+ only want one job per thread, so the number of jobs and the number of
+ threads are the same. */
+ gal_threads_spin_off(warp_wcsalign_init_convert, wa, wa->numthreads,
+ wa->numthreads, input->minmapsize,
+ input->quietmmap);
+
+ /* Pickle variables so other functions can access them */
+ wa->gcrn=gcrn;
+ wa->input=input;
+ wa->ncrn=4*es+4;
+ wa->v0=WARP_WCSALIGN_V0(es, os0, os1);
+
+ /* Determine the output image rotation direction so we can sort the
+ indices in counter clockwise order. This is necessary for the
+ 'gal_polygon_clip' function to work */
+ warp_check_output_orientation(wa);
+}
+
+
+
+
+
+void
+gal_warp_wcsalign_onpix(gal_warp_wcsalign_t *wa, size_t ind)
+{
+ size_t ic, temp, numinput=0;
+ gal_data_t *input=wa->input;
+ double xmin, xmax, ymin, ymax;
+ long xstart, ystart, xend, yend, x, y; /* Might be negative */
+ double filledarea, v, *ocrn=NULL, pcrn[8], opixarea;
+
+ size_t ncrn=wa->ncrn;
+ size_t is0=input->dsize[0];
+ size_t is1=input->dsize[1];
+ double *inputarr=input->array;
+ double *outputarr=wa->output->array;
+
+ size_t numcrn=0;
+ double ccrn[GAL_POLYGON_MAX_CORNERS], area;
+
+ /* Initialize the output pixel value: */
+ outputarr[ind] = filledarea = 0.0f;
+
+ if( wa->isccw==1 )
+ ocrn=warp_pixel_perimeter_cw(wa, ind);
+ else if( wa->isccw==0 )
+ ocrn=warp_pixel_perimeter_ccw(wa, ind);
+ else
+ error(EXIT_FAILURE, 0, "a bug! the code %d is not recognized as "
+ "a valid rotation orientation in "
+ "'gal_polygon_is_counterclockwise', this is not your fault, "
+ "something in the programming has gone wrong. Please contact "
+ "us at %s so we can correct it", wa->isccw, PACKAGE_BUGREPORT);
+
+ /* Find overlapping pixels */
+ xmin = DBL_MAX; ymin = DBL_MAX;
+ xmax = -DBL_MAX; ymax = -DBL_MAX;
+ for(ic=ncrn; ic--;)
+ {
+ temp=ic*2;
+ if(xmin > ocrn[ temp ]) { xmin = ocrn[ temp ]; }
+ if(xmax < ocrn[ temp ]) { xmax = ocrn[ temp ]; }
+ if(ymin > ocrn[ temp+1 ]) { ymin = ocrn[ temp+1 ]; }
+ if(ymax < ocrn[ temp+1 ]) { ymax = ocrn[ temp+1 ]; }
+ }
+
+ /* Start and end in both dimensions. */
+ xstart = GAL_DIMENSION_NEARESTINT_HALFHIGHER( xmin );
+ ystart = GAL_DIMENSION_NEARESTINT_HALFHIGHER( ymin );
+ xend = GAL_DIMENSION_NEARESTINT_HALFLOWER( xmax ) + 1;
+ yend = GAL_DIMENSION_NEARESTINT_HALFLOWER( ymax ) + 1;
+
+ /* Check which input pixels we are covering */
+ for(y=ystart;y<yend;++y)
+ {
+ /* If the pixel isn't in the image (note that the pixel
+ coordinates start from 1), skip this pixel. */
+ if( y<1 || y>is0 ) continue;
+
+ /* Y of base pixel vertices, in pixel coords. */
+ pcrn[1]=y-0.5f; pcrn[3]=y-0.5f;
+ pcrn[5]=y+0.5f; pcrn[7]=y+0.5f;
+
+ for(x=xstart;x<xend;++x)
+ {
+ if( x<1 || x>is1 ) continue;
+
+ /* X of base pixel vertices, in pixel coords. */
+ pcrn[0]=x-0.5f; pcrn[2]=x+0.5f;
+ pcrn[4]=x+0.5f; pcrn[6]=x-0.5f;
+
+ /* Read the value of the input pixel. */
+ v=inputarr[(y-1)*is1+x-1];
+
+ /* Find the overlapping (clipped) polygon: */
+ {
+ numcrn=0; /* initialize it. */
+ gal_polygon_clip(ocrn, ncrn, pcrn, 4, ccrn, &numcrn);
+ area=gal_polygon_area(ccrn, numcrn);
+
+ /* Add1 the fractional value of this pixel. If this output
+ pixel covers a NaN pixel in the input grid, then
+ calculate the area of this NaN pixel to account for it
+ later. */
+ if( !isnan(v) )
+ {
+ numinput+=1;
+ filledarea+=area;
+ outputarr[ind]+=v*area;
+
+ /* Check
+ printf("Check: numinput %zu filledarea %f "
+ "outputarr[%zu]=%f\n",
+ filledarea, ind, outputarr[ind]);
+ */
+ }
+ }
+ }
+ }
+
+ /* See if the pixel value should be set to NaN or not (because of not
+ enough coverage). Note that 'ocrn' is sorted in anti-clockwise
+ order already. */
+ opixarea=gal_polygon_area(ocrn, ncrn);
+ if( numinput && filledarea/opixarea < wa->coveredfrac-1e-5)
+ numinput=0;
+
+ /* Write the final value and return. */
+ if( numinput ==0) outputarr[ind]=NAN;
+
+ /* Clean up. */
+ free(ocrn);
+}
+
+
+
+
+
+void *
+gal_warp_wcsalign_onthread(void *inparam)
+{
+ size_t i, ind;
+ struct gal_threads_params *tprm=(struct gal_threads_params *)inparam;
+ gal_warp_wcsalign_t *wa=(gal_warp_wcsalign_t *)tprm->params;
+
+ /* Loop over pixels given from the 'warp' function */
+ for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
+ {
+ ind=tprm->indexs[i];
+ gal_warp_wcsalign_onpix(wa, ind);
+ }
+
+ /* Wait for all the other threads to finish, then return. */
+ if(tprm->b) { pthread_barrier_wait(tprm->b); }
+ return NULL;
+}
+
+
+
+
+
+/* Clean up the internally allocated variables from the 'wa' struct. */
+void
+gal_warp_wcsalign_free(gal_warp_wcsalign_t *wa)
+{
+ gal_list_data_free(wa->vertices);
+ wa->vertices=NULL;
+}
+
+
+
+
+
+/* Finalize the output 'gal_data_t' image in 'nl->output' */
+void
+gal_warp_wcsalign(gal_warp_wcsalign_t *nl)
+{
+ /* Calculate and allocate the output image size and WCS */
+ gal_warp_wcsalign_init(nl);
+
+ /* Fill the output image */
+ gal_threads_spin_off(gal_warp_wcsalign_onthread, nl, nl->output->size,
+ nl->numthreads, nl->input->minmapsize,
+ nl->input->quietmmap);
+
+ /* Clean up the internally allocated variables */
+ gal_warp_wcsalign_free(nl);
+}
diff --git a/lib/wcs.c b/lib/wcs.c
index ba83d6b5..1f98553a 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -483,6 +483,28 @@ gal_wcs_create(double *crpix, double *crval, double *cdelt,
+void
+gal_wcs_free(struct wcsprm *wcs)
+{
+ int status;
+
+ /* If it is already NULL, then don't continue. */
+ if(wcs==NULL) return;
+
+ /* Free the internal structure. */
+ status=wcsfree(wcs);
+ if(status)
+ error(EXIT_FAILURE, 0, "%s: WCSLIB wcsfree ERROR %d: %s",
+ __func__, status, wcs_errmsg[status]);
+
+ /* Free the actual structure. */
+ free(wcs);
+}
+
+
+
+
+
/* Extract the dimension name from CTYPE. */
char *
gal_wcs_dimension_name(struct wcsprm *wcs, size_t dimension)
@@ -1144,23 +1166,25 @@ gal_wcs_distortion_identify(struct wcsprm *wcs)
if( dispre != NULL )
{
- if( !strcmp(*dispre->dtype, "SIP") ) return GAL_WCS_DISTORTION_SIP;
- else if( !strcmp(*dispre->dtype, "TPD") ) return GAL_WCS_DISTORTION_TPD;
+ if( !strcmp(*dispre->dtype, "SIP") )
+ return GAL_WCS_DISTORTION_SIP;
+ else if( !strcmp(*dispre->dtype, "TPD") )
+ return GAL_WCS_DISTORTION_TPD;
else
- error(EXIT_FAILURE, 0, "%s: distortion '%s' isn't recognized in "
- "the 'dispre' structure of the given 'wcsprm'", __func__,
- *dispre->dtype);
+ return GAL_WCS_DISTORTION_INVALID;
}
else if( disseq != NULL )
{
- if( !strcmp(*disseq->dtype, "TPV") ) return GAL_WCS_DISTORTION_TPV;
- else if( !strcmp(*disseq->dtype, "TPD") ) return GAL_WCS_DISTORTION_TPD;
- else if( !strcmp(*disseq->dtype, "DSS") ) return GAL_WCS_DISTORTION_DSS;
- else if( !strcmp(*disseq->dtype, "WAT") ) return GAL_WCS_DISTORTION_WAT;
+ if( !strcmp(*disseq->dtype, "TPV") )
+ return GAL_WCS_DISTORTION_TPV;
+ else if( !strcmp(*disseq->dtype, "TPD") )
+ return GAL_WCS_DISTORTION_TPD;
+ else if( !strcmp(*disseq->dtype, "DSS") )
+ return GAL_WCS_DISTORTION_DSS;
+ else if( !strcmp(*disseq->dtype, "WAT") )
+ return GAL_WCS_DISTORTION_WAT;
else
- error(EXIT_FAILURE, 0, "%s: distortion '%s' isn't recognized in "
- "the 'disseq' structure of the given 'wcsprm'", __func__,
- *dispre->dtype);
+ return GAL_WCS_DISTORTION_INVALID;
}
/* Control should not reach here. */
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index b7c9459d..7b80f889 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -88,8 +88,6 @@ options=
-
-
# RUN THE PROCEDURES
# ==================
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 70693618: Warp: aligning the pixel and WCS coords, while correcting distortions,
Mohammad Akhlaghi <=