[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 58ec3af: Unified options used for single crop
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 58ec3af: Unified options used for single crop from center |
|
Date: |
Tue, 27 Jun 2017 18:15:27 -0400 (EDT) |
branch: master
commit 58ec3afdd2421b312ca689c1b709b4e3a16adc72
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Unified options used for single crop from center
Until now, if you wanted to crop a region around the pixel (A,B) in an
image with Crop, you had to use `--xc=A --yc=B'. On the other hand, if you
wanted to a crop around the RA and Dec of (R,D), you have to use `--ra=R
--dec=D'. Each standard also has its own separate option to specify the
width: `--iwidth' and `--wwidth'. In parallel, Crop also has the `--mode'
option to specify the coordinate standard (for cases like `--polygon',
where only one option is used for both standards).
This can all be confusing, just adds to the number of options (makes things
harder for the user) and also makes maintaining Crop harder. The main
reason for this is historical: Crop is Gnuastro's first program, so I
wasn't too experience and all these options grew very gradually.
This commit addresses these issues as described below. Since these changes
were significant, they have also been included in the NEWS file so users
downloading the source before the 0.4 version can be aware of these major
changes in an old program.
The new `--width' option is now used to define the width of a single
crop. Hence the old `--iwidth', `--wwidth' were removed. The units to
interpret the value to the option are specified by the `--mode'
option. With the new `--width' option it is also possible to define a
non-square crop (different widths along each dimension). In WCS mode, its
units are no longer arcseconds but are the same units of the WCS (degrees
for angles). `--width' can also accept fractions. So to set a width of 5
arcseconds, you can give it a value of `5/3600' for the angular
dimensions.
The new `--center' option is now used to define the center of a single
crop. Hence the old `--ra', `--dec', `--xc', `--yc' have been
removed. This new option can take multiple values (one value for each
dimension). Fractions are also acceptable.
The new `--coordcol' option is now used to determine the catalog columns
that define coordinates. Hence the old `--racol', `--deccol', `--xcol',
and `--ycol' have been removed. This new option can be called multiple
times and the order of its calling will be used for the column containing
the center in the respective dimension (in FITS format).
The following library functions have also been corrected.
`gal_box_overlap' now works on data of any dimensionality and thus also
needs the number of dimensions (elements in each input array).
`gal_box_border_from_center' now accepts an array of coordinates as one
argument and the number of dimensions as another. This allows it to work
on any dimensionality.
`gal_wcs_pixel_scale' now replaces the old `gal_wcs_pixel_scale_deg',
since it doesn't only apply to degrees. The pixel scale units are defined
by the units of the WCS.
---
NEWS | 51 ++
THANKS | 1 +
bin/crop/args.h | 174 ++----
bin/crop/astcrop.conf | 2 -
bin/crop/crop.c | 30 +-
bin/crop/main.h | 28 +-
bin/crop/onecrop.c | 248 ++++-----
bin/crop/onecrop.h | 11 +-
bin/crop/ui.c | 591 ++++++++++++++-------
bin/crop/ui.h | 20 +-
bin/crop/wcsmode.c | 273 +++++-----
bin/crop/wcsmode.h | 6 +-
bin/mkprof/mkprof.c | 7 +-
bin/warp/ui.c | 2 +-
bin/warp/warp.c | 2 +-
configure.ac | 2 +-
doc/gnuastro.texi | 314 +++++------
lib/blank.c | 2 +-
lib/box.c | 185 ++++---
lib/gnuastro-internal/options.h | 3 +
lib/gnuastro/box.h | 4 +-
lib/gnuastro/dimension.h | 4 +-
lib/gnuastro/wcs.h | 2 +-
lib/options.c | 22 +-
lib/wcs.c | 8 +-
tests/Makefile.am | 14 +-
tests/crop/imgcat.sh | 5 +-
tests/crop/{xcyc.sh => imgcenter.sh} | 5 +-
tests/crop/{xcycnoblank.sh => imgcenternoblank.sh} | 7 +-
tests/crop/section.sh | 2 +-
tests/crop/wcscat.sh | 5 +-
tests/crop/{radec.sh => wcscenter.sh} | 6 +-
tests/during-dev.sh | 8 +-
33 files changed, 1110 insertions(+), 934 deletions(-)
diff --git a/NEWS b/NEWS
index f231f95..b249c68 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,56 @@
GNU Astronomy Utilities NEWS -*- outline -*-
+* Noteworthy changes in release 0.3.XXX (library 2.0.0) (YYYY-MM-DD) [stable]
+
+** New features
+
+ `.fit' suffix is now also recognized as a FITS file extension in all
+ programs and libraries.
+
+ Arithmetic: now has a new `--globalhdu' (`-g') option which can be used
+ once for all the input images.
+
+** Removed features
+
+** Changed features
+
+ Crop: The new `--center' option is now used to define the center of a
+ single crop. Hence the old `--ra', `--dec', `--xc', `--yc' have been
+ removed. This new option can take multiple values (one value for each
+ dimension). Fractions are also acceptable.
+
+ Crop: The new `--width' option is now used to define the width of a
+ single crop. Hence the old `--iwidth', `--wwidth' were removed. The units
+ to interpret the value to the option are specified by the `--mode'
+ option. With the new `--width' option it is also possible to define a
+ non-square crop (different widths along each dimension). In WCS mode, its
+ units are no longer arcseconds but are the same units of the WCS (degrees
+ for angles). `--width' can also accept fractions. So to set a width of 5
+ arcseconds, you can give it a value of `5/3600' for the angular
+ dimensions.
+
+ Crop: The new `--coordcol' option is now used to determine the catalog
+ columns that define coordinates. Hence the old `--racol', `--deccol',
+ `--xcol', and `--ycol' have been removed. This new option can be called
+ multiple times and the order of its calling will be used for the column
+ containing the center in the respective dimension (in FITS format).
+
+ `gal_fits_img_info' now also returns the name and units of the dataset
+ (if they aren't NULL). So it takes two extra arguments.
+
+ `gal_box_overlap' now works on data of any dimensionality and thus also
+ needs the number of dimensions (elements in each input array).
+
+ `gal_box_border_from_center' now accepts an array of coordinates as one
+ argument and the number of dimensions as another. This allows it to work
+ on any dimensionality.
+
+ `gal_wcs_pixel_scale' now replaces the old `gal_wcs_pixel_scale_deg',
+ since it doesn't only apply to degrees. The pixel scale units are defined
+ by the units of the WCS.
+
+** Bug fixes
+
* Noteworthy changes in release 0.3 (library 1.0.0) (2017-06-01) [stable]
This is a full re-write of Gnuastro. Most importantly, Gnuastro now has a
diff --git a/THANKS b/THANKS
index 183b1d0..a8d34d3 100644
--- a/THANKS
+++ b/THANKS
@@ -31,6 +31,7 @@ support in Gnuastro. The list is ordered alphabetically.
Yahya Sefidbakht address@hidden
Richard Stallman address@hidden
Ole Streicher address@hidden
+ Ignacio Trujillo address@hidden
David Valls-Gabaud address@hidden
Christopher Willmer address@hidden
diff --git a/bin/crop/args.h b/bin/crop/args.h
index f973bd4..0cc9fbf 100644
--- a/bin/crop/args.h
+++ b/bin/crop/args.h
@@ -32,6 +32,20 @@ struct argp_option program_options[] =
{
/* Input. */
{
+ "mode",
+ UI_KEY_MODE,
+ "STR",
+ 0,
+ "Coordinate mode `img' or `wcs'",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->mode,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_parse_coordinate_mode
+ },
+ {
"hstartwcs",
UI_KEY_HSTARTWCS,
"INT",
@@ -108,7 +122,7 @@ struct argp_option program_options[] =
{
0, 0, 0, 0,
- "Crop by center (general settings)",
+ "Crop by center",
ARGS_GROUP_CENTER_GENERAL
},
{
@@ -125,92 +139,32 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "iwidth",
- UI_KEY_IWIDTH,
- "INT",
+ "width",
+ UI_KEY_WIDTH,
+ "FLT[,...]",
0,
- "Width (pixels) when crop defined by X,Y.",
+ "Width when crop is defined by its center.",
ARGS_GROUP_CENTER_GENERAL,
- &p->iwidthin,
- GAL_TYPE_SIZE_T,
- GAL_OPTIONS_RANGE_GT_0_ODD,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "wwidth",
- UI_KEY_WWIDTH,
- "FLT",
- 0,
- "Width (arcseconds) for crops defined by RA,Dec.",
- ARGS_GROUP_CENTER_GENERAL,
- &p->wwidth,
- GAL_TYPE_FLOAT64,
- GAL_OPTIONS_RANGE_GT_0,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
-
-
-
-
-
- {
- 0, 0, 0, 0,
- "Crop by center (single crop)",
- ARGS_GROUP_CENTER_SINGLE
- },
- {
- "ra",
- UI_KEY_RA,
- "FLT",
- 0,
- "Right ascension of one crop box center.",
- ARGS_GROUP_CENTER_SINGLE,
- &p->ra,
- GAL_TYPE_FLOAT64,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "dec",
- UI_KEY_DEC,
- "FLT",
- 0,
- "Declination of one crop box center.",
- ARGS_GROUP_CENTER_SINGLE,
- &p->dec,
- GAL_TYPE_FLOAT64,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "xc",
- UI_KEY_XC,
- "FLT",
- 0,
- "First axis position of one crop box center.",
- ARGS_GROUP_CENTER_SINGLE,
- &p->xc,
- GAL_TYPE_FLOAT64,
+ &p->width,
+ GAL_TYPE_STRING,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
+ GAL_OPTIONS_NOT_SET,
+ ui_parse_width_and_center
},
{
- "yc",
- UI_KEY_YC,
- "FLT",
+ "center",
+ UI_KEY_CENTER,
+ "FLT[,...]",
0,
- "Second axis position of one crop box center.",
- ARGS_GROUP_CENTER_SINGLE,
- &p->yc,
- GAL_TYPE_FLOAT64,
+ "Central coordinates of a single crop.",
+ ARGS_GROUP_CENTER_GENERAL,
+ &p->center,
+ GAL_TYPE_STRING,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
+ GAL_OPTIONS_NOT_SET,
+ ui_parse_width_and_center
},
@@ -221,7 +175,7 @@ struct argp_option program_options[] =
{
0, 0, 0, 0,
- "Crop by center (catalog)",
+ "Crop by center (when a catalog is given)",
ARGS_GROUP_CENTER_CATALOG
},
{
@@ -264,53 +218,14 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "racol",
- UI_KEY_RACOL,
- "STR/INT",
- 0,
- "Column number/info of Right Ascension (RA).",
- ARGS_GROUP_CENTER_CATALOG,
- &p->racol,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "deccol",
- UI_KEY_DECCOL,
- "STR/INT",
- 0,
- "Column number/info of Declination.",
- ARGS_GROUP_CENTER_CATALOG,
- &p->deccol,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "xcol",
- UI_KEY_XCOL,
+ "coordcol",
+ UI_KEY_COORDCOL,
"STR/INT",
0,
- "Column number/info of X (first FITS axis).",
+ "Column no./info containing coordinates.",
ARGS_GROUP_CENTER_CATALOG,
- &p->xcol,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "ycol",
- UI_KEY_YCOL,
- "STR/INT",
- 0,
- "Column number/info of Y (second FITS axis).",
- ARGS_GROUP_CENTER_CATALOG,
- &p->ycol,
- GAL_TYPE_STRING,
+ &p->coordcol,
+ GAL_TYPE_STRLL,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
@@ -371,19 +286,6 @@ struct argp_option program_options[] =
/* Operating mode */
- {
- "mode",
- UI_KEY_MODE,
- "STR",
- 0,
- "Coordinate mode `img' or `wcs'",
- GAL_OPTIONS_GROUP_OPERATING_MODE,
- &p->modestr,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
diff --git a/bin/crop/astcrop.conf b/bin/crop/astcrop.conf
index d58836a..b616160 100644
--- a/bin/crop/astcrop.conf
+++ b/bin/crop/astcrop.conf
@@ -19,8 +19,6 @@
# Input image and catalog parameters:
cathdu 1
- iwidth 201
- wwidth 3
hstartwcs 0
hendwcs 0
diff --git a/bin/crop/crop.c b/bin/crop/crop.c
index 0d12f95..db1dc3a 100644
--- a/bin/crop/crop.c
+++ b/bin/crop/crop.c
@@ -52,7 +52,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
truncated. In the log-file there is no truncation, therefore the log
file should be used for checking the outputs, not the outputs printed on
the screen. */
-void
+static void
crop_verbose_info(struct onecropparams *crp)
{
char *filestatus, *msg;
@@ -84,7 +84,7 @@ crop_verbose_info(struct onecropparams *crp)
/* Print final statistics in verbose mode. */
-void
+static void
crop_verbose_final(struct cropparams *p)
{
char *msg;
@@ -149,7 +149,7 @@ crop_verbose_final(struct cropparams *p)
-void
+static void
crop_write_to_log(struct onecropparams *crp)
{
char **strarr;
@@ -185,8 +185,8 @@ crop_write_to_log(struct onecropparams *crp)
-void *
-imgmodecrop(void *inparam)
+static void *
+crop_mode_img(void *inparam)
{
struct onecropparams *crp=(struct onecropparams *)inparam;
struct cropparams *p=crp->p;
@@ -210,7 +210,7 @@ imgmodecrop(void *inparam)
crp->out_ind=crp->indexs[i];
crp->outfits=NULL;
crp->numimg=1; /* In Image mode there is only one input image. */
- cropname(crp);
+ onecrop_name(crp);
/* Crop the image. */
onecrop(crp);
@@ -219,7 +219,7 @@ imgmodecrop(void *inparam)
if(crp->numimg)
{
/* Check if the center of the crop is filled or not. */
- crp->centerfilled=iscenterfilled(crp);
+ crp->centerfilled=onecrop_center_filled(crp);
/* Add the final headers and close output FITS image: */
gal_fits_key_write_version(crp->outfits, NULL, PROGRAM_STRING);
@@ -261,8 +261,8 @@ imgmodecrop(void *inparam)
-void *
-wcsmodecrop(void *inparam)
+static void *
+crop_mode_wcs(void *inparam)
{
struct onecropparams *crp=(struct onecropparams *)inparam;
struct cropparams *p=crp->p;
@@ -282,21 +282,21 @@ wcsmodecrop(void *inparam)
/* Set the sides of the crop in RA and Dec */
- setcsides(crp);
+ wcsmode_crop_corners(crp);
/* Go over all the images to see if this target is within their
range or not. */
crp->in_ind=0;
do
- if(radecoverlap(crp))
+ if(wcsmode_overlap(crp))
{
/* Open the input FITS file. */
crp->infits=gal_fits_hdu_open_format(p->imgs[crp->in_ind].name,
p->cp.hdu, 0);
/* If a name isn't set yet, set it. */
- if(crp->name==NULL) cropname(crp);
+ if(crp->name==NULL) onecrop_name(crp);
/* Increment the number of images used (necessary for the
header keywords that are written in `onecrop'). Then do the
@@ -322,7 +322,7 @@ wcsmodecrop(void *inparam)
if(crp->numimg)
{
/* See if the center is filled. */
- crp->centerfilled=iscenterfilled(crp);
+ crp->centerfilled=onecrop_center_filled(crp);
/* Write all the dependency versions and close the file. */
gal_fits_key_write_version(crp->outfits, NULL, PROGRAM_STRING);
@@ -340,7 +340,7 @@ wcsmodecrop(void *inparam)
}
else
{
- cropname(crp);
+ onecrop_name(crp);
crp->centerfilled=0;
}
@@ -399,7 +399,7 @@ crop(struct cropparams *p)
/* Set the function to run: */
- modefunction = p->mode==IMGCROP_MODE_IMG ? &imgmodecrop : &wcsmodecrop;
+ modefunction = p->mode==IMGCROP_MODE_IMG ? &crop_mode_img : &crop_mode_wcs;
/* Allocate the array of structures to keep the thread and parameters for
diff --git a/bin/crop/main.h b/bin/crop/main.h
index 48fd6ed..d9e3f88 100644
--- a/bin/crop/main.h
+++ b/bin/crop/main.h
@@ -40,7 +40,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Macros */
#define LOGFILENAME PROGRAM_EXEC".log"
#define FILENAME_BUFFER_IN_VERB 30
-
+#define MAXDIM 2
/* Modes of operation. */
@@ -56,8 +56,7 @@ enum crop_modes
/* The sides of the image keep the celestial coordinates of the four
- sides of this image. With respect to the pixels they are.
-*/
+ sides of this image. With respect to the pixels they are. */
struct inputimgs
{
char *name; /* File name of input image. */
@@ -81,43 +80,34 @@ struct cropparams
{
/* Directly from command-line */
struct gal_options_common_params cp; /* Common parameters. */
- struct gal_list_str_t *inputs; /* All input FITS files. */
+ gal_list_str_t *inputs; /* All input FITS files. */
size_t hstartwcs; /* Header keyword No. to start read WCS. */
size_t hendwcs; /* Header keyword No. to end read WCS. */
+ int mode; /* Image or WCS mode. */
uint8_t zeroisnotblank; /* ==1: In float or double, keep 0.0. */
uint8_t noblank; /* ==1: no blank (out of image) pixels. */
char *suffix; /* Ending of output file name. */
size_t checkcenter; /* width of a box to check for zeros */
- size_t iwidthin; /* Image mode width (in pixels). */
- double wwidth; /* WCS mode width (in arcseconds). */
- double ra; /* RA of one crop box center. */
- double dec; /* Dec of one crop box center. */
- double xc; /* Center point, one crop (FITS stnrd). */
- double yc; /* Center point, one crop (FITS stnrd). */
+ gal_data_t *center; /* Center position of crop. */
+ gal_data_t *width; /* Width of crop when defined by center. */
char *catname; /* Name of input catalog. */
char *cathdu; /* HDU of catalog if its a FITS file. */
char *namecol; /* Filename (without suffix) of crop col.*/
- char *racol; /* Catalog RA column */
- char *deccol; /* Catalog Dec column */
- char *xcol; /* Catalog X column */
- char *ycol; /* Catalog Y column */
+ gal_list_str_t *coordcol; /* Column in catalog with coordinates. */
char *section; /* Section string. */
char *polygon; /* Input string of polygon vertices. */
uint8_t outpolygon; /* ==1: Keep the inner polygon region. */
- char *modestr; /* ==1: will use X and Y coordiates. */
/* Internal */
- int mode; /* Image or WCS mode. */
size_t numin; /* Number of input images. */
size_t numout; /* Number of output images. */
- double *c1; /* First coordinate from catalog. */
- double *c2; /* Second coordinate from catalog. */
+ double **centercoords; /* The center coordinates. */
char **name; /* filename of crop in row. */
double *wpolygon; /* Array of WCS polygon vertices. */
double *ipolygon; /* Array of image polygon vertices. */
size_t nvertices; /* Number of polygon vertices. */
long iwidth[2]; /* Image mode width (in pixels). */
- double res; /* Resolution in arcseconds */
+ double *pixscale; /* Resolution in each dimension. */
time_t rawtime; /* Starting time of the program. */
int outnameisfile; /* Output filename is a directory. */
int type; /* Type of output(s). */
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index a4634f9..225de23 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -72,24 +72,30 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Read the section string and set the starting and ending pixels
based on that. */
void
-sectionparser(struct cropparams *p, size_t *dsize,
- long *fpixel, long *lpixel)
+onecrop_parse_section(struct cropparams *p, size_t *dsize,
+ long *fpixel, long *lpixel)
{
int add;
long read;
- size_t dim=0;
char *tailptr;
char forl='f', *pt=p->section;
long naxes[2]={dsize[1], dsize[0]};
+ size_t i, dim=0, ndim=p->imgs->ndim;
- /* If control reached here, then the cropped region is not defined by its
- center. So it makes no sense to check if the center is blank. */
+ /* When the user asks for a section of the dataset, then the cropped
+ region is not defined by its center. So it makes no sense to later
+ check if the center is blank or not. Hence, we will over-write it with
+ zero. */
p->checkcenter=0;
- /* Initialize the fpixel and lpixel arrays: */
- lpixel[0]=naxes[0];
- lpixel[1]=naxes[1];
- fpixel[0]=fpixel[1]=1;
+ /* Initialize the fpixel and lpixel arrays (note that `section' is only
+ defined in image mode, so there will only be one element in `imgs'. */
+ for(i=0;i<ndim;++i)
+ {
+ fpixel[i] = 1;
+ lpixel[i] = naxes[i] = p->imgs->dsize[ ndim - i - 1 ];
+ }
+
/* Parse the string. */
while(*pt!='\0')
@@ -99,7 +105,7 @@ sectionparser(struct cropparams *p, size_t *dsize,
{
case ',':
++dim;
- if(dim==2)
+ if(dim>=ndim)
error(EXIT_FAILURE, 0, "Extra `,` in `%s`", p->section);
forl='f';
++pt;
@@ -138,34 +144,29 @@ sectionparser(struct cropparams *p, size_t *dsize,
else continue;
}
- /* Put it in the correct array. Note that for the last point, we
- don't want to include the given pixel. Unlike CFITSIO, in Crop,
- the given intervals are not inclusive. fpixel and lpixel will be
- directly passed to CFITSIO. So we have to correct his here.*/
+ /* Put it in the correct array. */
if(forl=='f')
fpixel[dim] = add ? naxes[dim]+read : read;
else
lpixel[dim] = add ? naxes[dim]+read : read;
pt=tailptr;
-
- /* For a check:
- printf("\n\n[%ld, %ld]: fpixel=(%ld, %ld), lpixel=(%ld, %ld)\n\n",
- naxes[0], naxes[1],
- fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
- */
}
/* Make sure the first pixel is located before/below the last pixel. */
- if(fpixel[0]>lpixel[0] || fpixel[1]>lpixel[1])
- error(EXIT_FAILURE, 0, "the bottom left corner coordinates "
- "cannot be larger or equal to the top right's! Your section "
- "string (%s) has been read as: bottom left coordinate "
- "(%ld, %ld) to top right coordinate (%ld, %ld)",
- p->section, fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
-
- /*
- printf("\n%s\nfpixel=(%ld, %ld), lpixel=(%ld, %ld)\n\n", section,
- fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
+ for(i=0;i<ndim;++i)
+ if(fpixel[i]>lpixel[i])
+ error(EXIT_FAILURE, 0, "the bottom left corner coordinates "
+ "cannot be larger or equal to the top right's! Your section "
+ "string (%s) has been read as: bottom left coordinate "
+ "(%ld, %ld) to top right coordinate (%ld, %ld)",
+ p->section, fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
+
+ /* For a check:
+ printf("\n%s\n", p->section);
+ printf("fpixel: ("); for(i=0;i<ndim;++i) printf("%ld, ", fpixel[i]);
+ printf("\b\b)\n");
+ printf("lpixel: ("); for(i=0;i<ndim;++i) printf("%ld, ", lpixel[i]);
+ printf("\b\b)\n\n");
exit(0);
*/
}
@@ -176,7 +177,7 @@ sectionparser(struct cropparams *p, size_t *dsize,
void
-crop_polygonparser(struct cropparams *p)
+onecrop_parse_polygon(struct cropparams *p)
{
size_t dim=0;
char *tailptr;
@@ -281,8 +282,8 @@ crop_polygonparser(struct cropparams *p)
void
-imgpolygonflpixel(double *ipolygon, size_t nvertices, long *fpixel,
- long *lpixel)
+onecrop_ipolygon_fl(double *ipolygon, size_t nvertices, long *fpixel,
+ long *lpixel)
{
size_t i;
double minx=FLT_MAX, miny=FLT_MAX;
@@ -399,8 +400,8 @@ polygonmask(struct onecropparams *crp, void *array, long
*fpixel_i,
/*******************************************************************/
/****************** One crop. *********************/
/*******************************************************************/
-void
-changezerotonan(void *array, size_t size, int type)
+static void
+onecrop_zero_to_nan(void *array, size_t size, int type)
{
float *fp, *ffp;
double *dp, *fdp;
@@ -434,7 +435,7 @@ changezerotonan(void *array, size_t size, int type)
/* Set the output name and image output sizes. */
void
-cropname(struct onecropparams *crp)
+onecrop_name(struct onecropparams *crp)
{
char **strarr;
struct cropparams *p=crp->p;
@@ -477,73 +478,68 @@ cropname(struct onecropparams *crp)
-/* Find the first and last pixel of a crop from its center point (in
- image mode or WCS mode). */
-void
-cropflpixel(struct onecropparams *crp)
+/* Find the first and last pixel of a crop. */
+static void
+onecrop_flpixel(struct onecropparams *crp)
{
struct cropparams *p=crp->p;
- int ncoord=1, nelem=2, status[2]={0,0};
- size_t *dsize=p->imgs[crp->in_ind].dsize;
- double pixcrd[2], imgcrd[2], phi[1], theta[1];
+ size_t ndim=p->imgs->ndim;
+
+ double center[MAXDIM];
+ int ncoord=1, status=0;
+ size_t i, *dsize=p->imgs[crp->in_ind].dsize;
long *fpixel=crp->fpixel, *lpixel=crp->lpixel;
+ double pixcrd[MAXDIM], imgcrd[MAXDIM], phi[1], theta[1];
switch(p->mode)
{
case IMGCROP_MODE_IMG:
- if(p->catname)
- gal_box_border_from_center(p->c1[crp->out_ind], p->c2[crp->out_ind],
- p->iwidth, fpixel, lpixel);
- else if(!isnan(p->xc))
- gal_box_border_from_center(p->xc, p->yc, p->iwidth, fpixel, lpixel);
- else if(p->section)
- sectionparser(p, dsize, fpixel, lpixel);
- else if(p->polygon)
+ if(p->section) /* Defined by section. */
+ onecrop_parse_section(p, dsize, fpixel, lpixel);
+ else if(p->polygon) /* Defined by polygon. */
{
if(p->outpolygon==0)
- imgpolygonflpixel(p->ipolygon, p->nvertices, fpixel, lpixel);
+ onecrop_ipolygon_fl(p->ipolygon, p->nvertices, fpixel, lpixel);
}
else
- error(EXIT_FAILURE, 0, "a bug! In image mode, neither of the "
- "following has been set: a catalog, a central pixel, "
- "a section or a polygon in the image. Please contact us "
- "to see how it got to this impossible place! You should "
- "have been warned of this condition long before Crop "
- "reaches this point");
+ {
+ for(i=0;i<ndim;++i) center[i] = p->centercoords[i][crp->out_ind];
+ gal_box_border_from_center(center, ndim, p->iwidth, fpixel, lpixel);
+ }
break;
+
case IMGCROP_MODE_WCS: /* In wcsmode, crp->world is already filled. */
if(p->polygon) /* Note: p->iwidth was set based on p->wwidth. */
{
/* Fill crp->ipolygon in wcspolygonpixel, then set flpixel*/
fillcrpipolygon(crp);
if(p->outpolygon==0)
- imgpolygonflpixel(crp->ipolygon, p->nvertices, fpixel, lpixel);
+ onecrop_ipolygon_fl(crp->ipolygon, p->nvertices, fpixel, lpixel);
}
else
{
- if(wcss2p(p->imgs[crp->in_ind].wcs, ncoord, nelem, crp->world,
- phi, theta, imgcrd, pixcrd, status) )
- if(status[0] || status[1])
+ /* Convert `crp->world' (in WCS) into `pixcrd' (image coord). */
+ if(wcss2p(p->imgs[crp->in_ind].wcs, ncoord, ndim, crp->world,
+ phi, theta, imgcrd, pixcrd, &status) )
+ if(status)
error(EXIT_FAILURE, 0, "%s: wcss2p error %d: %s", __func__,
- status[0] ? status[0] : status[1],
- wcs_errmsg[status[0] ? status[0] : status[1]]);
- gal_box_border_from_center(pixcrd[0], pixcrd[1], p->iwidth, fpixel,
- lpixel);
- /*
- printf("\n(%f, %f): (%ld, %ld) -- (%ld, %ld)\n\n", pixcrd[0],
- pixcrd[1], fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
- */
+ status, wcs_errmsg[status]);
+
+ /* Find the first and last pixels of this crop. */
+ gal_box_border_from_center(pixcrd, ndim, p->iwidth, fpixel, lpixel);
}
break;
+
default:
error(EXIT_FAILURE, 0, "%s: a bug! The domain (WCS or image) are not "
"set. Please contact us at %s so we can see how it got to this "
"impossible place", __func__, PACKAGE_BUGREPORT);
}
+
/* If the user only wants regions outside to the polygon, then set
the fpixel and lpixel to cover the full input image. */
if(p->polygon && p->outpolygon)
@@ -568,17 +564,17 @@ cropflpixel(struct onecropparams *crp)
blank areas or not. While the fpixel_i and lpixel_i arrays keep the
first and last pixels after the blank pixels have been removed.
*/
-void
-firstcropmakearray(struct onecropparams *crp, long *fpixel_i,
+static void
+onecrop_make_array(struct onecropparams *crp, long *fpixel_i,
long *lpixel_i, long *fpixel_c, long *lpixel_c)
{
- size_t i;
+ double crpix;
fitsfile *ofp;
- long naxes[2];
- int type=crp->p->type;
- double crpix0, crpix1;
- int naxis=2, status=0;
+ long naxes[MAXDIM];
char *outname=crp->name;
+ char cpname[FLEN_KEYWORD];
+ int status, type=crp->p->type;
+ size_t i, ndim=crp->p->imgs->ndim;
char *cp, *cpf, blankrec[80], titlerec[80];
char startblank[]=" / ";
struct inputimgs *img=&crp->p->imgs[crp->in_ind];
@@ -593,16 +589,14 @@ firstcropmakearray(struct onecropparams *crp, long
*fpixel_i,
/* Set the size of the output, in WCS mode, noblank==0. */
if(crp->p->noblank && crp->p->mode==IMGCROP_MODE_IMG)
- {
- fpixel_c[0]=fpixel_c[1]=1;
- lpixel_c[0]=naxes[0]=lpixel_i[0]-fpixel_i[0]+1;
- lpixel_c[1]=naxes[1]=lpixel_i[1]-fpixel_i[1]+1;
- }
+ for(i=0;i<ndim;++i)
+ {
+ fpixel_c[i] = 1;
+ lpixel_c[i] = naxes[i]=lpixel_i[i]-fpixel_i[i]+1;
+ }
else
- {
- naxes[0]=crp->lpixel[0]-crp->fpixel[0]+1;
- naxes[1]=crp->lpixel[1]-crp->fpixel[1]+1;
- }
+ for(i=0;i<ndim;++i)
+ naxes[i]=crp->lpixel[i] - crp->fpixel[i] + 1;
/* Create the FITS file with a blank first extension, then close it, so
@@ -615,13 +609,15 @@ firstcropmakearray(struct onecropparams *crp, long
*fpixel_i,
fits_create_img(ofp, SHORT_IMG, 0, naxes, &status);
fits_close_file(ofp, &status);
+
/* Create the output crop image. */
fits_open_file(&crp->outfits, outname, READWRITE, &status);
fits_create_img(crp->outfits, gal_fits_type_to_bitpix(type),
- naxis, naxes, &status);
+ ndim, naxes, &status);
gal_fits_io_error(status, "creating image");
ofp=crp->outfits;
+
/* When CFITSIO creates a FITS extension it adds two comments linking to
the FITS paper. Since we are mentioning the version of CFITSIO and
only use its ruitines to read/write from/to FITS files, this is
@@ -633,7 +629,8 @@ firstcropmakearray(struct onecropparams *crp, long
*fpixel_i,
fits_delete_key(ofp, "COMMENT", &status);
status=0;
- /* Write the blank value if necessary. */
+
+ /* Write the blank value as a FITS keyword if necessary. */
if( type!=GAL_TYPE_FLOAT32 && type!=GAL_TYPE_FLOAT64 )
if(fits_write_key(ofp, gal_fits_type_to_datatype(crp->p->type), "BLANK",
crp->p->bitnul, "pixels with no data", &status) )
@@ -646,8 +643,7 @@ firstcropmakearray(struct onecropparams *crp, long
*fpixel_i,
update the header keywords. */
if(img->wcs)
{
- crpix0 = img->wcs->crpix[0] - (fpixel_i[0]-1) + (fpixel_c[0]-1);
- crpix1 = img->wcs->crpix[1] - (fpixel_i[1]-1) + (fpixel_c[1]-1);
+ /* Write the WCS title and common WCS information. */
if(fits_write_record(ofp, blankrec, &status))
gal_fits_io_error(status, NULL);
sprintf(titlerec, "%sWCS information", startblank);
@@ -656,9 +652,17 @@ firstcropmakearray(struct onecropparams *crp, long
*fpixel_i,
fits_write_record(ofp, titlerec, &status);
for(i=0;i<img->nwcskeys-1;++i)
fits_write_record(ofp, &img->wcstxt[i*80], &status);
- fits_update_key(ofp, TDOUBLE, "CRPIX1", &crpix0, NULL, &status);
- fits_update_key(ofp, TDOUBLE, "CRPIX2", &crpix1, NULL, &status);
gal_fits_io_error(status, NULL);
+
+ /* Correct the CRPIX keywords. */
+ for(i=0;i<ndim;++i)
+ {
+ sprintf(cpname, "CRPIX%zu", i+1);
+ crpix = img->wcs->crpix[i] - (fpixel_i[i]-1) + (fpixel_c[i]-1);
+ fits_update_key(ofp, TDOUBLE, cpname, &crpix, NULL, &status);
+ gal_fits_io_error(status, NULL);
+ }
+
}
@@ -677,12 +681,7 @@ firstcropmakearray(struct onecropparams *crp, long
*fpixel_i,
/* The starting and ending points are set in the onecropparams structure
- for one crop from one image. Crop that region out.
-
- return values are:
- 0: No crop was made (not in range of image).
- 1: The input image covered at least part of the crop image.
- */
+ for one crop from one image. Crop that region out of the input. */
void
onecrop(struct onecropparams *crp)
{
@@ -690,39 +689,45 @@ onecrop(struct onecropparams *crp)
struct inputimgs *img=&p->imgs[crp->in_ind];
void *array;
- size_t cropsize;
int status=0, anynul=0;
char basename[FLEN_KEYWORD];
fitsfile *ifp=crp->infits, *ofp;
gal_fits_list_key_t *headers=NULL;
- long fpixel_o[2], lpixel_o[2], inc[2]={1,1};
+ size_t i, j, cropsize=1, ndim=img->ndim;
char region[FLEN_VALUE], regionkey[FLEN_KEYWORD];
- long naxes[]={img->dsize[1], img->dsize[0]}, fpixel_i[2] , lpixel_i[2];
+ long fpixel_o[MAXDIM], lpixel_o[MAXDIM], inc[MAXDIM];
+ long naxes[MAXDIM], fpixel_i[MAXDIM] , lpixel_i[MAXDIM];
+
+ /* Fill the `naxes' and `inc' arrays. */
+ for(i=0;i<ndim;++i)
+ {
+ inc[ i ] = 1;
+ naxes[ i ] = img->dsize[ ndim - i - 1 ];
+ }
/* Find the first and last pixel of this crop box from this input
- image. If the outer polygon region is to be kept, then set the
- sides to the image sides.*/
- cropflpixel(crp);
- fpixel_i[0]=crp->fpixel[0]; fpixel_i[1]=crp->fpixel[1];
- lpixel_i[0]=crp->lpixel[0]; lpixel_i[1]=crp->lpixel[1];
+ image. Then copy the first and last pixels into the `_i' arrays.*/
+ onecrop_flpixel(crp);
+ memcpy(fpixel_i, crp->fpixel, ndim*sizeof *fpixel_i);
+ memcpy(lpixel_i, crp->lpixel, ndim*sizeof *lpixel_i);
+
/* Find the overlap and apply it if there is any overlap. */
- if( gal_box_overlap(naxes, fpixel_i, lpixel_i, fpixel_o, lpixel_o) )
+ if( gal_box_overlap(naxes, fpixel_i, lpixel_i, fpixel_o, lpixel_o, ndim) )
{
/* Make the output FITS image and initialize it with an array of
- NaN or BLANK values. Note that for FLOAT_IMG and DOUBLE_IMG,
- it will automatically fill them with the NaN value.*/
+ NaN or BLANK values. */
if(crp->outfits==NULL)
- firstcropmakearray(crp, fpixel_i, lpixel_i, fpixel_o, lpixel_o);
+ onecrop_make_array(crp, fpixel_i, lpixel_i, fpixel_o, lpixel_o);
ofp=crp->outfits;
- /* Read the desired part of the image, then write it into this
- array. */
+ /* Allocate an array to keep the desired crop region, then read
+ the desired pixels onto it. */
status=0;
- cropsize=(lpixel_i[0]-fpixel_i[0]+1)*(lpixel_i[1]-fpixel_i[1]+1);
+ for(i=0;i<ndim;++i) cropsize *= ( lpixel_i[i] - fpixel_i[i] + 1 );
array=gal_data_malloc_array(p->type, cropsize, __func__, "array");
if(fits_read_subset(ifp, gal_fits_type_to_datatype(p->type),
fpixel_i, lpixel_i, inc, p->bitnul, array,
@@ -736,7 +741,7 @@ onecrop(struct onecropparams *crp)
if(p->zeroisnotblank==0
&& (p->type==GAL_TYPE_FLOAT32
|| p->type==GAL_TYPE_FLOAT64) )
- changezerotonan(array, cropsize, p->type);
+ onecrop_zero_to_nan(array, cropsize, p->type);
/* If a polygon is given, remove all the pixels within or
@@ -759,14 +764,19 @@ onecrop(struct onecropparams *crp)
gal_fits_io_error(status, NULL);
+ /* Write the selected region of this image as a string to include as
+ a FITS keyword. Then we want to delete the last coma `,'.*/
+ j=0;
+ for(i=0;i<ndim;++i)
+ j += sprintf(®ion[j], "%ld:%ld,", fpixel_i[i], lpixel_i[i]);
+ region[j-1]='\0';
+
+
/* A section has been added to the cropped image from this input
- image, so increment crp->imgcount and save the information of
- this image. */
+ image, so save the information of this image. */
sprintf(basename, "ICF%zu", crp->numimg);
gal_fits_key_write_filename(basename, img->name, &headers);
sprintf(regionkey, "%sPIX", basename);
- sprintf(region, "%ld:%ld,%ld:%ld", fpixel_i[0], lpixel_i[0],
- fpixel_i[1], lpixel_i[1]);
gal_fits_key_list_add_end(&headers, GAL_TYPE_STRING, regionkey,
0, region, 0, "Range of pixels used for "
"this output.", 0, NULL);
@@ -807,7 +817,7 @@ onecrop(struct onecropparams *crp)
/****************** Check center *********************/
/*******************************************************************/
int
-iscenterfilled(struct onecropparams *crp)
+onecrop_center_filled(struct onecropparams *crp)
{
struct cropparams *p=crp->p;
diff --git a/bin/crop/onecrop.h b/bin/crop/onecrop.h
index 2384567..59e5177 100644
--- a/bin/crop/onecrop.h
+++ b/bin/crop/onecrop.h
@@ -52,7 +52,7 @@ struct onecropparams
/* For log */
char *name; /* Filename of crop. */
size_t numimg; /* Number of images used to make this crop. */
- unsigned char centerfilled; /* ==1 if the center is filled. */
+ unsigned char centerfilled; /* ==1 if the center is filled. */
/* Thread parameters. */
size_t *indexs; /* Indexs to be used in this thread. */
@@ -60,19 +60,16 @@ struct onecropparams
};
void
-crop_polygonparser(struct cropparams *p);
+onecrop_parse_polygon(struct cropparams *p);
void
-cropname(struct onecropparams *crp);
-
-void
-cropflpixel(struct onecropparams *crp);
+onecrop_name(struct onecropparams *crp);
void
onecrop(struct onecropparams *crp);
int
-iscenterfilled(struct onecropparams *crp);
+onecrop_center_filled(struct onecropparams *crp);
void
crop_print_log(struct onecropparams *p);
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 9b8cf8e..a8cc827 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/blank.h>
#include <gnuastro/table.h>
+#include <gnuastro/dimension.h>
#include <gnuastro-internal/timing.h>
#include <gnuastro-internal/options.h>
@@ -64,7 +65,7 @@ const char *
argp_program_bug_address = PACKAGE_BUGREPORT;
static char
-args_doc[] = "[ASCIIcatalog] ASTRdata ...";
+args_doc[] = "[Crop-identifiers] ASTRdata ...";
const char
doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will create cutouts, "
@@ -133,7 +134,6 @@ ui_initialize_options(struct cropparams *p,
/* Initalize necessary parameters. */
- p->xc=p->yc=p->ra=p->dec=NAN;
p->mode = IMGCROP_MODE_INVALID;
cp->searchin = GAL_TABLE_SEARCH_INVALID;
@@ -217,14 +217,106 @@ parse_opt(int key, char *arg, struct argp_state *state)
+/* Parse the mode to interpret the given coordinates. */
+void *
+ui_parse_coordinate_mode(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *junk)
+{
+ char *outstr;
+ /* We want to print the stored values. */
+ if(lineno==-1)
+ {
+ gal_checkset_allocate_copy( *(int *)(option->value)==IMGCROP_MODE_IMG
+ ? "img" : "wcs", &outstr );
+ return outstr;
+ }
+ else
+ {
+ if (!strcmp(arg, "img")) *(int *)(option->value)=IMGCROP_MODE_IMG;
+ else if (!strcmp(arg, "wcs")) *(int *)(option->value)=IMGCROP_MODE_WCS;
+ else
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "`%s' (value to "
+ "`--mode') not recognized as an input mode. "
+ "Recognized values are `img' and `wcs'. This option "
+ "is necessary to identify the nature of your input "
+ "coordinates.\n\n"
+ "Please run the following command for more "
+ "information (press the `SPACE' key to go down and "
+ "`q' to return to the command-line):\n\n"
+ " $ info gnuastro \"Crop modes\"\n", arg);
+ return NULL;
+ }
+}
+/* Parse the width and center coordinates from the comman-line or
+ configuration files. */
+void *
+ui_parse_width_and_center(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *junk)
+{
+ size_t i, nc;
+ double *darray;
+ gal_data_t *values;
+ char *str, sstr[GAL_OPTIONS_STATIC_MEM_FOR_VALUES];
+ /* We want to print the stored values. */
+ if(lineno==-1)
+ {
+ /* Set the value pointer to the correct type. */
+ values = *(gal_data_t **)(option->value);
+ darray = values->array;
+ /* Write the values into a string. */
+ nc=0;
+ for(i=0;i<values->size;++i)
+ {
+ if( nc > GAL_OPTIONS_STATIC_MEM_FOR_VALUES-100 )
+ error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so we "
+ "can address the problem. The number of necessary "
+ "characters in the statically allocated string has become "
+ "too close to %d", __func__, PACKAGE_BUGREPORT,
+ GAL_OPTIONS_STATIC_MEM_FOR_VALUES);
+ nc += sprintf(sstr+nc, "%g,", darray[i]);
+ }
+ sstr[nc-1]='\0';
+
+ /* Copy the string into a dynamically allocated space, because it
+ will be freed later.*/
+ gal_checkset_allocate_copy(sstr, &str);
+ return str;
+ }
+ else
+ {
+ /* If the option is already set, then ignore it. */
+ if(option->set) return NULL;
+
+ /* Read the values. */
+ values=gal_options_parse_list_of_numbers(arg, filename, lineno);
+ *(gal_data_t **)(option->value) = values;
+
+ /* If we are on the width option, then make sure none of the values
+ are negative or zero. */
+ if(!strcmp(option->name, "width"))
+ {
+ darray=values->array;
+ for(i=0;i<values->size;++i)
+ if(darray[i]<=0.0f)
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "%g is <=0. "
+ "The values to the `--width' option must be "
+ "larger than zero. The complete input to this "
+ "option was `%s' (%g is input number %zu)",
+ darray[i], arg, darray[i], i+1);
+ }
+
+ /* The return value is only for printing mode. */
+ return NULL;
+ }
+}
@@ -233,68 +325,40 @@ parse_opt(int key, char *arg, struct argp_state *state)
-/**************************************************************/
-/*************** Sanity Check *******************/
-/**************************************************************/
-/* Read and check ONLY the options. When arguments are involved, do the
- check in `ui_check_options_and_arguments'. */
-static void
-ui_read_check_only_options(struct cropparams *p)
-{
- int checksum;
- /* Read the mode from the string the user specified. */
- if(p->modestr)
- {
- if (!strcmp(p->modestr, "img")) p->mode=IMGCROP_MODE_IMG;
- else if (!strcmp(p->modestr, "wcs")) p->mode=IMGCROP_MODE_WCS;
- else
- error(EXIT_FAILURE, 0, "`%s' (value to `--mode') not recognized as "
- "an input mode. Recognized values are `img' and `wcs'. This "
- "option is necessary when the inputs are not sufficient to "
- "identify the nature of the coordinates.\n\n"
- "Please run the following command for more information "
- "(press the `SPACE' key to go down and `q' to return to the "
- "command-line):\n\n"
- " $ info gnuastro \"Crop modes\"\n", p->modestr);
- }
- /* Check that if one of the coordinates is given, the other is also
- given. Note that if both are not given, their sum will be 0, if both
- are given, then the sum will be 2. If only one is given, then the sum
- will be 1. */
- if( (isnan(p->ra) + isnan(p->dec)) == 1 )
- error(EXIT_FAILURE, 0, "no `--%s' given, it is mandatory with the `--%s' "
- "option", isnan(p->ra)?"ra":"dec", isnan(p->ra)?"dec":"ra");
- if( (isnan(p->xc) + isnan(p->yc)) == 1 )
- error(EXIT_FAILURE, 0, "no `--%s' given, it is mandatory with the `--%s' "
- "option", isnan(p->xc)?"xc":"yc", isnan(p->xc)?"yc":"xc");
- if( ((p->racol!=NULL) + (p->deccol!=NULL)) == 1 )
- error(EXIT_FAILURE, 0, "no `--%scol' given, it is mandatory with the "
- "`--%scol' option", p->racol?"dec":"ra",
- p->racol?"ra":"dec");
- if( ((p->xcol!=NULL) + (p->ycol!=NULL)) == 1 )
- error(EXIT_FAILURE, 0, "no `--%scol' given, it is mandatory with the "
- "`--%scol' option", p->xcol?"y":"x", p->xcol?"x":"y");
- /* Make sure that the single crop modes are not called together. */
- checksum = ( (!isnan(p->xc)) + (!isnan(p->ra)) + (p->catname!=NULL)
- + (p->section!=NULL) + (p->polygon!=NULL) );
+
+/**************************************************************/
+/*************** Sanity Check *******************/
+/**************************************************************/
+/* Read and check ONLY the options. When arguments are involved, do the
+ check in `ui_check_options_and_arguments'. */
+static void
+ui_read_check_only_options(struct cropparams *p)
+{
+ int checksum;
+
+ /* Make sure that only one of the crop definitions is given. */
+ checksum = ( (p->center!=NULL)
+ + (p->catname!=NULL)
+ + (p->section!=NULL)
+ + (p->polygon!=NULL) );
switch(checksum)
{
case 0:
error(EXIT_FAILURE, 0, "no crop definition. You can use any of the "
- "following options to define the crop(s): (`--xc', `--yc'), "
- "(`--ra', `--dec'), `--catalog', `--section', `--polygon'. "
- "Please run this command for more information:\n\n"
+ "following options to define the crop(s): `--center', "
+ "`--catalog', `--section', or `--polygon'. Please run this "
+ "command for more information:\n\n"
" $ info gnuastro \"Crop modes\"\n");
case 1:
/* Everything is ok, just ignore the switch structure. */
@@ -302,16 +366,21 @@ ui_read_check_only_options(struct cropparams *p)
default:
error(EXIT_FAILURE, 0, "more than one crop type specified. In each "
"run, only one crop definition is acceptable on the "
- "command-line, or in configuration files. You have called "
- "%s%s%s%s%s\b\b.",
- !isnan(p->xc) ? "(`--xc', `--yc'), " : "",
- !isnan(p->ra) ? "(`--ra', `--dec'), " : "",
+ "command-line or in configuration files. You have called: "
+ "%s%s%s%s\b\b.",
+ p->center!=NULL ? "`--center', " : "",
p->catname!=NULL ? "`--catalog', " : "",
p->section!=NULL ? "`--section', " : "",
p->polygon!=NULL ? "`--polygon', " : "");
}
+ /* Section is currentlyl only defined in Image mode. */
+ if(p->section && p->mode!=IMGCROP_MODE_IMG)
+ error(EXIT_FAILURE, 0, "The `--section' option is only available in "
+ "image coordinate mode, currently it doesn't work with WCS mode. "
+ "Please run with `--mode=img' and change the values accordingly");
+
/* Sanity checks and mode setting based on the desired crop. */
if(p->catname)
@@ -334,50 +403,22 @@ ui_read_check_only_options(struct cropparams *p)
/* Atleast one of the (X,Y), and (RA,Dec) set of columns are
necessary. Note that we have checked that they are together if
given, so we only need to check one of the two in each couple. */
- if(p->xcol==NULL && p->racol==NULL)
- error(EXIT_FAILURE, 0, "no crop center position column given "
- "to read from the input catalog (`%s'). Please use either of "
- "(`--xcol', `--ycol') or (`--racol', `--deccol'). For more "
- "information on how to select columns in Gnuastro, please run "
- "the following command:\n\n"
+ if(p->coordcol==NULL)
+ error(EXIT_FAILURE, 0, "no crop center columns given to read from "
+ "the input catalog (`%s'). Please use `--coordcol' two times "
+ "to specify the column keeping the center position the "
+ "respective dimension.\n\n"
+ "For more information on how to select columns in Gnuastro, "
+ "please run the following command:\n\n"
" $ info gnuastro \"Selecting table columns\"", p->catname);
-
- /* If only image columns are specified, then we have Image mode, if
- only WCS columns are specified, we have WSC mode. If both are
- specified, the mode is mandatory */
- if(p->xcol && p->racol)
- {
- if(p->mode==IMGCROP_MODE_INVALID)
- error(EXIT_FAILURE, 0, "both image and WCS coordinate columns "
- "are specified to read the center of the crops in the "
- "input catalog (%s). You can use the `--mode=img', or "
- "`--mode=wcs' options to specify which set of columns "
- "should be used", p->catname);
- }
- else
- p->mode = p->xcol ? IMGCROP_MODE_IMG : IMGCROP_MODE_WCS;
- }
- else if(p->polygon)
- {
- if(p->mode==IMGCROP_MODE_INVALID)
- error(EXIT_FAILURE, 0, "the polygon option works in both image and "
- "wcs mode, please specify how the vertices should be "
- "interpreted with `--mode=img', or `--mode=wcs' options to "
- "specify which set of columns should be used");
}
- else if(!isnan(p->xc))
- p->mode = IMGCROP_MODE_IMG;
- else if(!isnan(p->ra))
- p->mode = IMGCROP_MODE_WCS;
- else if( p->section)
- p->mode = IMGCROP_MODE_IMG;
/* Parse the polygon vertices if they are given to make sure that it is
in the proper format. */
if(p->polygon)
{
- crop_polygonparser(p);
+ onecrop_parse_polygon(p);
if(p->nvertices<3)
error(EXIT_FAILURE, 0, "a polygon has to have 3 or more vertices, "
"you have only given %zu (%s)", p->nvertices, p->polygon);
@@ -405,12 +446,11 @@ ui_read_check_only_options(struct cropparams *p)
static void
ui_check_options_and_arguments(struct cropparams *p)
{
- /* Make sure we do actually have inputs. */
+ /* Make sure we actually have inputs. */
if(p->inputs==NULL)
error(EXIT_FAILURE, 0, "no input file given");
- /* Make sure an input file name was given and if it was a FITS file, that
- a HDU is also given. */
+ /* Make sure that a HDU is also given. */
if(p->cp.hdu==NULL )
error(EXIT_FAILURE, 0, "no HDU specified. When the input is a FITS "
"file, a HDU must also be specified, you can use the `--hdu' "
@@ -475,88 +515,186 @@ ui_check_options_and_arguments(struct cropparams *p)
/**************************************************************/
/*************** Preparations *******************/
/**************************************************************/
+/* When the crop is defined by its center, the final width that we need
+ must be in actual number of pixels (an integer). But the user's values
+ can be in WCS mode or even in image mode, they may be non-integers. */
+static void
+ui_set_iwidth(struct cropparams *p)
+{
+ gal_data_t *newwidth;
+ double pwidth, *warray;
+ size_t i, ndim=p->imgs->ndim;
+
+
+ /* Make sure a width value is actually given. */
+ if(p->width==NULL)
+ error(EXIT_FAILURE, 0, "no crop width specified. When crops are "
+ "defined by their center (with `--center' or `--catalog') a "
+ "width is necessary (using the `--width' option)");
+
+ /* Make sure that the width array only has one element or the same number
+ of elements as the input's dimensions. */
+ if(p->width->size!=ndim && p->width->size!=1)
+ error(EXIT_FAILURE, 0, "%zu values give to `--width', but input is %zu "
+ "dimensions. It can only take either one value (same width in all "
+ "dimensions), or the same number as the input's dimensions",
+ p->width->size, ndim);
+
+ /* If the width array has only one value, that single value should be
+ used for all dimensions. */
+ if(p->width->size==1)
+ {
+ /* Allocate the new width dataset. */
+ newwidth=gal_data_alloc(NULL, p->width->type, 1, &ndim, NULL, 0, -1,
+ NULL, NULL, NULL);
+
+ /* Fill the new width. */
+ warray=newwidth->array;
+ for(i=0;i<ndim;++i) warray[i] = *(double *)(p->width->array);
+
+ /* Free the old (single element) width dataset and put the new one in
+ its place. */
+ gal_data_free(p->width);
+ p->width=newwidth;
+ }
+ else warray=p->width->array;
+
+ /* Fill in `p->iwidth' depending on the mode. */
+ for(i=0;i<ndim;++i)
+ {
+ /* Set iwidth. */
+ if(p->mode==IMGCROP_MODE_WCS)
+ {
+ /* Convert the width in units of the input's WCS into pixels. */
+ pwidth = warray[i]/p->pixscale[i];
+ if(pwidth<3 || pwidth>50000)
+ error(EXIT_FAILURE, 0, "%g (width along dimension %zu) "
+ "translates to %.0f pixels. This is probably not what "
+ "you wanted. Note that the resolution in this dimension "
+ "is %g", warray[i], i+1, pwidth, p->pixscale[i]);
+
+ /* Write the single valued width in WCS and the image width for
+ this dimension. */
+ p->iwidth[i]=GAL_DIMENSION_FLT_TO_INT(pwidth);
+ if(p->iwidth[i]%2==0)
+ {
+ p->iwidth[i] += 1;
+ warray[i] += p->pixscale[i];
+ }
+ }
+ else
+ {
+ p->iwidth[i]=GAL_DIMENSION_FLT_TO_INT(warray[i]);
+ if(p->iwidth[i]%2==0) p->iwidth[i] += 1;
+ }
+ }
+
+ /* For a check:
+ printf("Width: "); for(i=0;i<ndim;++i) printf("\t%ld\n\t", p->iwidth[i]);
+ exit(0);
+ */
+}
+
+
+
+
+
static void
ui_read_cols(struct cropparams *p)
{
- char *colname;
- size_t counter=0;
- int toomanycols=0;
+ char colname[100];
gal_list_str_t *colstrs=NULL;
gal_data_t *cols, *tmp, *corrtype=NULL;
- char *ax1col = p->mode==IMGCROP_MODE_IMG ? p->xcol : p->racol;
- char *ax2col = p->mode==IMGCROP_MODE_IMG ? p->ycol : p->deccol;
+ size_t ncoordcols, counter=0, dcounter=0, ndim=p->imgs->ndim;
+
- /* Specify the order of columns. */
+ /* See if the number of columns given for coordinates corresponds to the
+ number of dimensions of the input dataset. */
+ if(p->coordcol)
+ {
+ /* Check if the number of columns given for coordinates is the same
+ as the number of dimensions in the input dataset(s). */
+ ncoordcols=gal_list_str_number(p->coordcol);
+ if( ncoordcols != ndim)
+ error(EXIT_FAILURE, 0, "`--coordcol' was called %zu times, but the "
+ "input dataset%s %zu dimensions, these values must not be "
+ "different. Recall that through `--coordcol' you are "
+ "specifying the columns containing the coordinates of the "
+ "center of the crop in a catalog", ncoordcols,
+ (p->numin==1?" has":"s have"), ndim);
+ }
+ else
+ error(EXIT_FAILURE, 0, "no coordinate columns specified. When a catalog"
+ "is given, it is necessary to identify which columns identify "
+ "the coordinate values in which dimension.\n\n"
+ "You can do this by calling `--coordcol' multiple times, the "
+ "order must be in the same order as the input's dimensions. "
+ "For more information on how to select columns in Gnuastro, "
+ "please run the following command:\n\n"
+ " $ info gnuastro \"Selecting table columns\"");
+
+
+ /* If a name column was also given, the put that as the first column to
+ read, otherwise just use the given set of columns (in the same
+ order). */
if(p->namecol)
- gal_list_str_add(&colstrs, p->namecol, 0);
- gal_list_str_add(&colstrs, ax1col, 0);
- gal_list_str_add(&colstrs, ax2col, 0);
+ {
+ gal_list_str_add(&colstrs, p->namecol, 0);
+ colstrs->next=p->coordcol;
+ }
+ else colstrs=p->coordcol;
+
/* Read the desired columns from the file. */
cols=gal_table_read(p->catname, p->cathdu, colstrs, p->cp.searchin,
p->cp.ignorecase, p->cp.minmapsize);
- /* Set the number of objects. */
+
+ /* Set the number of objects (rows in each column). */
p->numout=cols->size;
- /* For a sanity check, make sure that the total number of columns read is
- the same as those that were wanted (it might be more). */
+
+ /* Make sure more columns were not read (the name matchings might result
+ in more than one column being read from the inputs). */
+ if( gal_list_data_number(cols) != ndim + (p->namecol!=NULL) )
+ gal_tableintern_error_col_selection(p->catname, p->cathdu, "too many "
+ "columns were selected by the given "
+ "values to the options ending in "
+ "`col'.");
+
+
+ /* Put the information in each column in the proper place. */
while(cols!=NULL)
{
/* Pop out the top node. */
tmp=gal_list_data_pop(&cols);
- /* Note that the input was a linked list, so the output order is the
- inverse of the input order. For the position, we will store the
- values into the `x' and `y' arrays even if they are RA/Dec. */
+ /* See which column it is. */
switch(++counter)
{
- case 3:
- colname="crop name prefix";
+ case 1:
if(p->namecol)
{
- if(tmp->type==GAL_TYPE_STRING)
- {
- p->name=tmp->array;
- tmp->array=NULL;
- gal_data_free(tmp);
- }
- else
- {
- corrtype=gal_data_copy_to_new_type_free(tmp,
- GAL_TYPE_STRING);
- p->name=corrtype->array;
- }
+ sprintf(colname, "crop name prefix");
+ corrtype = (tmp->type==GAL_TYPE_STRING ? tmp
+ : gal_data_copy_to_new_type_free(tmp,
+ GAL_TYPE_STRING));
+ p->name=corrtype->array;
}
else
- toomanycols=1;
- break;
-
- case 2:
- colname="first axis position";
- corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
- p->c1=corrtype->array;
- break;
-
- case 1:
- colname="second axis position";
- corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
- p->c2=corrtype->array;
+ {
+ sprintf(colname, "position in dimension %zu", dcounter+1);
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+ p->centercoords[ dcounter++ ]=corrtype->array;
+ }
break;
- /* If the index isn't recognized, then it is larger, showing that
- there was more than one match for the given criteria */
default:
- toomanycols=1;
+ sprintf(colname, "position in dimension %zu", dcounter+1);
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+ p->centercoords[ dcounter++ ] = corrtype->array;
}
- /* Print an error if there were too many columns. */
- if(toomanycols)
- gal_tableintern_error_col_selection(p->catname, p->cathdu, "too many "
- "columns were selected by the "
- "given values to the options "
- "ending in `col'.");
-
/* Sanity check and clean up. Note that it might happen that the
input structure is already freed. In that case, `corrtype' will be
NULL. */
@@ -564,8 +702,10 @@ ui_read_cols(struct cropparams *p)
{
/* Make sure there are no blank values in this column. */
if( gal_blank_present(corrtype, 1) )
- error(EXIT_FAILURE, 0, "%s column has blank values. "
- "Input columns cannot contain blank values", colname);
+ error(EXIT_FAILURE, 0, "%s: column with %s has blank values. "
+ "Input columns must not contain blank values",
+ gal_fits_name_save_as_string(p->catname, p->cathdu),
+ colname);
/* Free the unnecessary sturcture information. The correct-type
(`corrtype') data structure's array is necessary for later
@@ -583,6 +723,47 @@ ui_read_cols(struct cropparams *p)
+static void
+ui_prepare_center(struct cropparams *p)
+{
+ double *carray;
+ size_t i, ndim=p->imgs->ndim;
+
+ /* Allocate space to keep the central positions. */
+ errno=0;
+ p->centercoords=malloc(ndim*sizeof *p->centercoords);
+ if( p->centercoords==NULL )
+ error(EXIT_FAILURE, 0, "%s: %zu bytes for `p->centercoords'",
+ __func__, ndim*sizeof *p->centercoords);
+
+
+ /* Set the integer widths of the crop(s) when defined by center. */
+ ui_set_iwidth(p);
+
+ /* For a catalog, we have a separate function, but for a single center
+ value, put the center values into an array. This will essentially
+ simulate a catalog with one row. So from this point on, there is no
+ difference between a catalog input and a central position input. */
+ if(p->catname)
+ ui_read_cols(p);
+ else
+ {
+ carray=p->center->array;
+ for(i=0;i<ndim;++i)
+ {
+ p->centercoords[i]=gal_data_malloc_array(GAL_TYPE_FLOAT64,
+ 1, __func__,
+ "p->centercoords[i]");
+ p->centercoords[i][0]=carray[i];
+ }
+ }
+}
+
+
+
+
+
+
/* Add all the columns of the log file. Just note that since this is a
linked list, we have to add them in the opposite order. */
static void
@@ -590,10 +771,10 @@ ui_make_log(struct cropparams *p)
{
char *comment;
- /* Return if no long file is to be created. */
+ /* Return if no long file was requested. */
if(p->cp.log==0) return;
- /* If central pixels are filled. */
+ /* Column to specify if the central pixels are filled. */
asprintf(&comment, "Are the central pixels filled? (1: yes, 0: no, "
"%u: not checked)", GAL_BLANK_UINT8);
gal_list_data_add_alloc(&p->log, NULL, GAL_TYPE_UINT8, 1, &p->numout,
@@ -601,12 +782,12 @@ ui_make_log(struct cropparams *p)
"bool", comment);
free(comment);
- /* Number of images used. */
+ /* Column for number of datasets used in this crop. */
gal_list_data_add_alloc(&p->log, NULL, GAL_TYPE_UINT16, 1, &p->numout,
NULL, 1, p->cp.minmapsize, "NUM_INPUTS", "count",
- "Number of input images used to make this crop.");
+ "Number of input datasets used to make this crop.");
- /* Row number in input catalog. */
+ /* Filename of crop. */
gal_list_data_add_alloc(&p->log, NULL, GAL_TYPE_STRING, 1, &p->numout,
NULL, 1, p->cp.minmapsize, "CROP_NAME", "name",
"File name of crop.");
@@ -619,16 +800,10 @@ ui_make_log(struct cropparams *p)
void
ui_preparations(struct cropparams *p)
{
- char *msg;
fitsfile *tmpfits;
- struct timeval t1;
- size_t input_counter;
struct inputimgs *img;
int status, firsttype=0;
-
-
- /* Set the initial iwidth. */
- p->iwidth[0]=p->iwidth[1]=p->iwidthin;
+ size_t input_counter, firstndim=0;
/* For polygon and section, there should be no center checking. */
@@ -636,22 +811,6 @@ ui_preparations(struct cropparams *p)
p->checkcenter=0;
- /* Read the columns if there is a catalog, otherwise, set the number
- of images (crops) to 1.*/
- if(p->catname)
- ui_read_cols(p);
- else
- p->numout=1;
-
-
- /* Everything is ready, notify the user of the program starting. */
- if(!p->cp.quiet)
- {
- gettimeofday(&t1, NULL);
- printf(PROGRAM_NAME" started on %s", ctime(&p->rawtime));
- }
-
-
/* Allocate space for all the input images. This is done here because
WCSLIB is unfortunately not thread-safe when reading the WCS
information from the FITS files. In cases where the number of cropped
@@ -663,8 +822,8 @@ ui_preparations(struct cropparams *p)
errno=0;
p->imgs=malloc(p->numin*sizeof *p->imgs);
if(p->imgs==NULL)
- error(EXIT_FAILURE, errno, "ui.c: %zu bytes for p->imgs",
- p->numin*sizeof *p->imgs);
+ error(EXIT_FAILURE, errno, "%s: %zu bytes for p->imgs",
+ __func__, p->numin*sizeof *p->imgs);
/* Fill in the WCS information of each image. */
@@ -691,23 +850,37 @@ ui_preparations(struct cropparams *p)
else
if(p->mode==IMGCROP_MODE_WCS)
error(EXIT_FAILURE, 0, "the WCS structure of %s (hdu: %s) "
- "image is not recognized. So RA and Dec cannot be used "
- "as input. You can try with pixel coordinates in the "
- "Image Mode (note that the crops will lack WCS "
- "header information)", img->name, p->cp.hdu);
+ "image is not recognized. So WCS mode cannot be used "
+ "as input coordinates. You can try with pixel coordinates "
+ "with `--mode=img'", img->name, p->cp.hdu);
fits_close_file(tmpfits, &status);
gal_fits_io_error(status, NULL);
- /* Make sure all the images have the same type. */
+ /* Make sure all the images have the same type and dimensions. */
if(firsttype==0)
{
- firsttype=p->type;
+ /* Set the basic information. */
+ firsttype = p->type;
+ firstndim = img->ndim;
p->bitnul = gal_blank_alloc_write(p->type);
+
+ /* Make sure the number of dimensions is supported. */
+ if(firstndim>MAXDIM)
+ error(EXIT_FAILURE, 0, "%s: is as %zu dimensional dataset, Crop "
+ "currently only supports a maximum of %d dimensions",
+ img->name, firstndim, MAXDIM);
+
+ /* Make sure the number of coordinates given for center
+ correspond to the dimensionality of the data. */
+ if(p->center && p->center->size!=firstndim)
+ error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, but "
+ "%zu coordinates were given to `--center'", img->name,
+ p->cp.hdu, firstndim, p->center->size);
}
else
{
if(firsttype!=p->type)
- error(EXIT_FAILURE, 0, "%s: type is `%s' while revious image(s) "
+ error(EXIT_FAILURE, 0, "%s: type is `%s' while previous input(s) "
"were `%s' type. All inputs must have the same pixel data "
"type.\n\nYou can use Gnuastro's Arithmetic program to "
"convert `%s' to `%s', please run this command for more "
@@ -717,21 +890,26 @@ ui_preparations(struct cropparams *p)
img->name, gal_type_name(p->type, 1),
gal_type_name(firsttype, 1), img->name,
gal_type_name(p->type, 1));
+ if(firstndim!=img->ndim)
+ error(EXIT_FAILURE, 0, "%s: type has %zu dimensions, while "
+ "previous input(s) had %zu dimensions. All inputs must "
+ "have the same number of dimensions", img->name, img->ndim,
+ firstndim);
}
- /* In WCS mode, Check resolution and get the first pixel
- positions. */
- if(p->mode==IMGCROP_MODE_WCS) wcs_check_prepare(p, img);
+ /* In WCS mode, we need some additional preparations. */
+ if(p->mode==IMGCROP_MODE_WCS) wcsmode_check_prepare(p, img);
}
- /* Report timing: */
- if(!p->cp.quiet)
- {
- asprintf(&msg, "Read metadata of %zu image%s.", p->numin,
- p->numin>1 ? "s" : "");
- gal_timing_report(&t1, msg, 1);
- }
+ /* Unify central crop methods into `p->centercoords'. */
+ if(p->catname || p->center)
+ ui_prepare_center(p);
+
+
+ /* `ui_read_cols' set the number of output crops when a catalog was
+ given, in all other cases, we only have one output. */
+ if(p->catname==NULL) p->numout=1;
/* Prepare the log file if the user has asked for it. */
@@ -756,6 +934,7 @@ ui_preparations(struct cropparams *p)
+
/**************************************************************/
/************ Set the parameters *************/
/**************************************************************/
@@ -763,6 +942,8 @@ ui_preparations(struct cropparams *p)
void
ui_read_check_inputs_setup(int argc, char *argv[], struct cropparams *p)
{
+ char *msg;
+ struct timeval t1;
struct gal_options_common_params *cp=&p->cp;
@@ -808,8 +989,29 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct
cropparams *p)
ui_check_options_and_arguments(p);
+ /* To see how long it takes to read meta-data. */
+ if(!p->cp.quiet) gettimeofday(&t1, NULL);
+
+
/* Read/allocate all the necessary starting arrays. */
ui_preparations(p);
+
+
+ /* Report timing: */
+ if(!p->cp.quiet)
+ {
+ printf(PROGRAM_NAME" started on %s", ctime(&p->rawtime));
+ asprintf(&msg, "Read metadata of %zu dataset%s.", p->numin,
+ p->numin>1 ? "s" : "");
+ gal_timing_report(&t1, msg, 1);
+ if(p->numout>1)
+ {
+ asprintf(&msg, "Will try making %zu crops (from catalog).",
+ p->numout);
+ gal_timing_report(NULL, msg, 1);
+ }
+ }
+
}
@@ -840,8 +1042,7 @@ ui_free_report(struct cropparams *p, struct timeval *t1)
size_t i;
/* Free the simple arrays (if they were set). */
- if(p->c1) free(p->c1);
- if(p->c2) free(p->c2);
+ gal_data_free(p->center);
if(p->cp.hdu) free(p->cp.hdu);
if(p->cathdu) free(p->cathdu);
if(p->catname) free(p->catname);
diff --git a/bin/crop/ui.h b/bin/crop/ui.h
index fdb9a70..46d520f 100644
--- a/bin/crop/ui.h
+++ b/bin/crop/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- e k m t u v
+ a d e f g i j k m r t u v y
A B E G H J L Q R W X Y
*/
enum option_keys_enum
@@ -37,23 +37,15 @@ enum option_keys_enum
/* With short-option version. */
UI_KEY_CATALOG = 'C',
UI_KEY_NOBLANK = 'b',
- UI_KEY_CHECKCENTER = 'c',
UI_KEY_SUFFIX = 'p',
UI_KEY_NAMECOL = 'n',
- UI_KEY_RACOL = 'f',
- UI_KEY_DECCOL = 'g',
- UI_KEY_RA = 'r',
- UI_KEY_DEC = 'd',
- UI_KEY_XCOL = 'i',
- UI_KEY_YCOL = 'j',
- UI_KEY_XC = 'x',
- UI_KEY_YC = 'y',
- UI_KEY_IWIDTH = 'a',
- UI_KEY_WWIDTH = 'w',
UI_KEY_SECTION = 's',
UI_KEY_POLYGON = 'l',
UI_KEY_ZEROISNOTBLANK = 'z',
UI_KEY_MODE = 'O',
+ UI_KEY_WIDTH = 'w',
+ UI_KEY_CENTER = 'c',
+ UI_KEY_COORDCOL = 'x',
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
@@ -61,12 +53,16 @@ enum option_keys_enum
UI_KEY_HSTARTWCS,
UI_KEY_HENDWCS,
UI_KEY_OUTPOLYGON,
+ UI_KEY_CHECKCENTER,
};
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct cropparams *p);
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index 34179a9..e443db4 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <error.h>
#include <stdio.h>
#include <float.h>
+#include <string.h>
#include <stdlib.h>
#include <gnuastro/wcs.h>
@@ -42,17 +43,20 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
+
+
/*******************************************************************/
/**************** Check for ui.c *********************/
/*******************************************************************/
/* This function is called from ui.c. Its job is to check the WCS
values of this */
void
-wcs_check_prepare(struct cropparams *p, struct inputimgs *img)
+wcsmode_check_prepare(struct cropparams *p, struct inputimgs *img)
{
- double twidth, *pixscale;
+ double *pixscale;
struct wcsprm *wcs=img->wcs;
- int i, status[4]={0,0,0,0}, ncoord=4, nelem=2;
+ int i, status[4]={0,0,0,0}, ncorners=4;
+ size_t *dsize=img->dsize, ndim=img->ndim;
double imgcrd[8], phi[4], theta[4], pixcrd[8];
@@ -81,93 +85,91 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs
*img)
img->name, p->cp.hdu);
- /* Check the image resolution (if the pixels are actually a square), then
- compare the resolution with the other input images. Due to floating
- point errors, some very small differences might exist in the pixel
- scale, so break out with an error only if the pixel scales are more
- different than 1e-6. */
- pixscale=gal_wcs_pixel_scale_deg(wcs);
+ /* Check the nature of the coordinates, currently we can only support RA
+ and Dec, other modes haven't been checked. */
+ if( strcmp(wcs->ctype[0], "RA---TAN")
+ || strcmp(wcs->ctype[1], "DEC--TAN") )
+ error(EXIT_FAILURE, 0, "currently the only WCS types usable are "
+ "`RA---TAN' and `DEC--TAN' for the first and second axises "
+ "respectively. The WCS types of `%s' (hdu %s) are `%s' and `%s' "
+ "respectively", img->name, p->cp.hdu, wcs->ctype[0], wcs->ctype[1]);
+
+
+ /* Check if the pixels are actually a square, then compare the resolution
+ with the other input images. Due to floating point errors, some very
+ small differences might exist in the pixel scale, so break out with an
+ error only if the pixel scales are more different than 1e-6. */
+ pixscale=gal_wcs_pixel_scale(wcs);
if( fabs(pixscale[0]-pixscale[1])/pixscale[0] > 1e-6 )
error(EXIT_FAILURE, 0, "%s: HDU %s: The pixel scale along "
"the two image axises is not the same. The first axis "
"is %.15g deg/pixel, while the second is %.15g",
img->name, p->cp.hdu, pixscale[1], pixscale[0]);
- if(p->res==0.0f)
+ if(p->pixscale)
{
- /* Set the resolution of the image. */
- p->res=pixscale[0];
+ for(i=0;i<ndim;++i)
+ if(p->pixscale[i] != pixscale[i])
+ error(EXIT_FAILURE, 0, "%s (hdu %s): has resolution of %g along "
+ "dimension %d. However, previously checked input(s) had "
+ "a resolution of %g in this dimension", img->name, p->cp.hdu,
+ pixscale[i], i+1, p->pixscale[i]);
free(pixscale);
-
- /* Set the widths such that iwidth and wwidth are exactly the same
- (within their different units ofcourse). Also make sure that the
- image size is an odd number (so the central pixel is in the
- center). */
- p->wwidth/=3600; /* Convert the width to degrees. */
- twidth=p->wwidth/p->res;
- if(twidth<3)
- error(EXIT_FAILURE, 0, "--wwidth = %f (arcseconds) translates "
- "to %.0f pixels in scale of input image(s). This is probably "
- "not what you want", p->wwidth*3600, twidth);
- p->iwidth[0] = (twidth-(long)twidth)>0.5 ? twidth+1 : twidth;
- if(p->iwidth[0]%2==0)
- {
- p->iwidth[0]+=1;
- p->wwidth+=p->res;
- }
- p->iwidth[1]=p->iwidth[0];
}
else
- {
- if(p->res!=pixscale[0])
- error(EXIT_FAILURE, 0, "%s: HDU %s: The resolution of "
- "this image is %f arcseconds/pixel while the "
- "previously checked input image(s) had a resolution "
- "of %f", img->name, p->cp.hdu, 3600*wcs->pc[3],
- 3600*p->res);
- free(pixscale);
- }
+ p->pixscale=pixscale;
+
+
+ /* Set the coordinates of the dataset's corners. Note that `dsize' is in
+ C order, while pixcrd is in FITS order.*/
+ pixcrd[0] = 1; pixcrd[1] = 1;
+ pixcrd[2] = dsize[1]; pixcrd[3] = 1;
+ pixcrd[4] = 1; pixcrd[5] = dsize[0];
+ pixcrd[6] = dsize[1]; pixcrd[7] = dsize[0];
- /* Get the coordinates of the first pixel in the image. Note that `dsize'
- is in C axises, while pixcrd is in FITS axises. */
- pixcrd[0]=1; pixcrd[1]=1;
- pixcrd[2]=img->dsize[1]; pixcrd[3]=1;
- pixcrd[4]=1; pixcrd[5]=img->dsize[0];
- pixcrd[6]=img->dsize[1]; pixcrd[7]=img->dsize[0];
- wcsp2s(wcs, ncoord, nelem, pixcrd, imgcrd, phi, theta,
+ /* Get the coordinates of the corners of the dataset in WCS. */
+ wcsp2s(wcs, ncorners, ndim, pixcrd, imgcrd, phi, theta,
img->corners, status);
+
/* Check if there was no error in the conversion. */
- for(i=0;i<4;++i)
+ for(i=0;i<ncorners;++i)
if(status[i])
error(EXIT_FAILURE, 0, "wcsp2s ERROR %d in row %d of pixcrd: %s",
i, status[i], wcs_errmsg[status[i]]);
- /* Fill in the size of the image in celestial degrees from the first
- pixel in the image. Note that `dsize' is in C axises, while pixcrd is
- in FITS axises. */
- img->sized[0]=img->dsize[1]*p->res/cos(img->corners[1]*M_PI/180);
- img->sized[1]=img->dsize[0]*p->res;
+ /* Fill in the size of the dataset in WCS from the first pixel in the
+ image. Note that `dsize' is in C axises, while the `pixscale',
+ `corners' and `sized' are in FITS axises. */
+ img->sized[0] = ( img->dsize[1] * p->pixscale[0]
+ / cos( img->corners[1] * M_PI / 180 ) );
+ img->sized[1] = img->dsize[0] * p->pixscale[1];
+
+ /* In case the image crosses the equator, we will calculate these values
+ here so later on, we don't have to calculate them on every check. See
+ the explanation above `point_in_dataset'.
- /* In case the image crosses the equator, we will calculate these
- values here so later on, we don't have to calculate them on every
- check. See the explanation above radecoverlap.*/
- if( img->corners[1]*(img->corners[1]+img->sized[1]) < 0 )
+ Note that in both 2D and 3D data, the declination is in the second
+ coordinate (index 1). */
+ if( img->corners[1] * (img->corners[1]+img->sized[1]) < 0 )
{
- /* re in the explanations. */
+ /* re in the comments above `point_in_dataset'. */
img->equatorcorr[0]=img->corners[0]
-0.5*img->sized[0]*(1-cos(img->corners[1]*M_PI/180));
- /* sre in the explanations. */
+ /* sre in the comments above `point_in_dataset'. */
img->equatorcorr[1]=img->sized[0]*cos(img->corners[1]*M_PI/180);
}
/* Just to check:
- printf("\n\n%s:\n(%.10f, %.10f)\n(%.10f, %.10f)"
- "\n(%.10f, %.10f)\n(%.10f, %.10f)\n\n", img->name,
+ printf("\n\n%s:\n", img->name);
+ printf("(%.10f, %.10f)\n"
+ "(%.10f, %.10f)\n"
+ "(%.10f, %.10f)\n"
+ "(%.10f, %.10f)\n\n",
img->corners[0], img->corners[1],
img->corners[2], img->corners[3],
img->corners[4], img->corners[5],
@@ -200,19 +202,21 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs
*img)
/*******************************************************************/
/* Set the four sides around the point of interest in RA and Dec.
- NOTE: In this format we are working on here (where the image is
- aligned with the celestial coordinates), the declination is
- measured on a great circle, while the right ascension is not. So we
- have to consider the
+ NOTE: When the image is aligned with the celestial coordinates (current
+ working paradigm), the declination is measured on a great circle, while
+ the right ascension is not. So we have to consider this in calculating
+ the difference in RA.
*/
void
-setcsides(struct onecropparams *crp)
+wcsmode_crop_corners(struct onecropparams *crp)
{
- size_t i;
- double h, hr, r, d, dr; /* The second r is for radians. */
struct cropparams *p=crp->p;
+
+ size_t i, ndim=p->imgs->ndim;
double minra=FLT_MAX, mindec=FLT_MAX;
double maxra=-FLT_MAX, maxdec=-FLT_MAX;
+ double r, d, dr, h[MAXDIM], hr[MAXDIM];
+ size_t rmini=-1, rmaxi=-1, dmini=-1, dmaxi=-1;
/* Set the four corners of the WCS region. */
if(p->polygon)
@@ -234,58 +238,75 @@ setcsides(struct onecropparams *crp)
}
else
{
- if(p->catname)
- {
- r=crp->world[0]=p->c1[crp->out_ind];
- d=crp->world[1]=p->c2[crp->out_ind];
- }
- else
- {
- r=crp->world[0]=p->ra;
- d=crp->world[1]=p->dec;
- }
+ /* Set the RA and Dec to use as center. */
+ r=crp->world[0]=p->centercoords[0][crp->out_ind];
+ d=crp->world[1]=p->centercoords[1][crp->out_ind];
+
- h=p->wwidth/2;
+ /* Calculate the declination in radians for easy readability. */
dr=d*M_PI/180;
- hr=h*M_PI/180;
- /* Set the four corners of this crop. */
- crp->corners[0] = r+h/cos(dr-hr); crp->corners[1] = d-h; /* Bt Lf */
- crp->corners[2] = r-h/cos(dr-hr); crp->corners[3] = d-h; /* Bt Rt */
- crp->corners[4] = r+h/cos(dr+hr); crp->corners[5] = d+h; /* Tp Lf */
- crp->corners[6] = r-h/cos(dr+hr); crp->corners[7] = d+h; /* Tp Rt */
+ /* Set the half width in each dimension. For the angular dimensions,
+ also calculate it in radians. */
+ hr[0] = ( h[0] = ((double *)(p->width->array))[0] / 2 ) * M_PI / 180;
+ hr[1] = ( h[1] = ((double *)(p->width->array))[1] / 2 ) * M_PI / 180;
+
+ /* Set the corners of this crop. */
+ crp->corners[0] = r+h[0]/cos(dr-hr[1]);
+ crp->corners[1] = d-h[1]; /* Bottom left. */
+
+ crp->corners[2] = r-h[0]/cos(dr-hr[1]);
+ crp->corners[3] = d-h[1]; /* Bottom Right. */
+
+ crp->corners[4] = r+h[0]/cos(dr+hr[1]);
+ crp->corners[5] = d+h[1]; /* Top Left. */
+
+ crp->corners[6] = r-h[0]/cos(dr+hr[1]);
+ crp->corners[7] = d+h[1]; /* Top Right. */
}
- /* Set the bottom width and height of the crop in degrees. Note that
- the width changes as the height changes, so here we want the
- height and the lowest declination. Note that on the bottom edge,
- corners[0] is the maximum RA and corners[2] is the minimum RA.
- For all the region, corners[5] is one of the maximum declinations
- and corners[3] is one of the the minimum declinations.*/
- crp->sized[0]=( (crp->corners[0]-crp->corners[2])
- / cos(crp->corners[1]*M_PI/180) );
- crp->sized[1]=crp->corners[5]-crp->corners[3];
+ /* Set the bottom width and height of the crop in degrees. Note that the
+ width changes as the height changes, so here we want the height and
+ the lowest declination. Note that in 2D on the bottom edge, corners[0]
+ is the maximum RA and corners[2] is the minimum RA. For all the 2D
+ region, corners[5] is one of the maximum declinations and corners[1]
+ is one of the the minimum declinations.
+
+ North and south hemispheres are no problem: When using the center,
+ they are set properly (in any hemisphere) and for a polygon, the
+ minimums and maximums are automatically found. */
+ rmini = ndim; /* First element in second corner. */
+ rmaxi = 0; /* First element. */
+ dmini = 1; /* Second element. */
+ dmaxi = 5; /* Second element in third corner. */
+ crp->sized[0]=( (crp->corners[rmaxi]-crp->corners[rmini])
+ / cos(crp->corners[dmini]*M_PI/180) );
+ crp->sized[1]=crp->corners[dmaxi]-crp->corners[dmini];
+
/* In case the crop crosses the equator, then we need these two
- corrections. See the complete explanations below. */
+ corrections. See the complete explanations above `point_in_dataset'. */
if(crp->corners[1]*(crp->corners[1]+crp->sized[1]) < 0 )
{
- /* re in the explanations. */
+ /* re in the explanations above `point_in_dataset'. */
crp->equatorcorr[0]=crp->corners[0]
-0.5*crp->sized[0]*(1-cos(crp->corners[1]*M_PI/180));
- /* sre in the explanations. */
+ /* sre in the explanations above `point_in_dataset'. */
crp->equatorcorr[1]=crp->sized[0]*cos(crp->corners[1]*M_PI/180);
}
/* Just to check:
- printf("\n\nCorner 1: (%.10f, %.10f)\n"
- "Corner 2: (%.10f, %.10f)\nCorner 3: (%.10f, %.10f)\n"
- "Corner 4: (%.10f, %.10f)\n\n",
+ printf("\n\n%g, %g:\n", r, d);
+ printf("\t(%.10f, %.10f)\n"
+ "\t(%.10f, %.10f)\n"
+ "\t(%.10f, %.10f)\n"
+ "\t(%.10f, %.10f)\n\n",
crp->corners[0], crp->corners[1],
crp->corners[2], crp->corners[3],
crp->corners[4], crp->corners[5],
crp->corners[6], crp->corners[7]);
+ exit(0);
*/
}
@@ -434,21 +455,29 @@ fillcrpipolygon(struct onecropparams *crp)
(South) equations of above as before and for those points that have
a positive declination, we use the North formula but replacing r1
with re, d1 with 0 and sr with sre.
-*/
-/*
- p[2]: RA and Dec of a point (rp and dp above).
- i[2]: RA and Dec of first point in image. (r1 and d1 above).
- s[2]: Size of the image in degrees (sr and sd above).
- c[2]: Corrections if equator is passed, (se and sre above).
-*/
-int
-radecinimg(double *p, double *i, double *s, double *c)
+
+
+ INPUTS
+ ------
+ p[]: Point coordinates (rp and dp above).
+ i[]: Coordinates of first pixel in image. (r1 and d1 above).
+ s[]: Size/width of box (sr and sd above).
+ c[]: Corrections if equator is passed, (se and sre above).
+
+
+ IMPORTANT: It is assumed that the dimensions are ordered with:
+
+ 0: RA
+ 1: Dec
+ 2: Third dimension (independent of RA and Dec). */
+static int
+point_in_dataset(double *p, double *i, double *s, double *c, size_t ndim)
{
double n;
- /* First check the declination. If it is not in range, you can
- safely return 0.*/
- if(p[1]>=i[1] && p[1]<=i[1]+s[1])
+ /* In the RA and Dec checks, first check the declination. If it is not in
+ range, you can safely return 0. */
+ if(p[1]>=i[1] && p[1]<=i[1]+s[1])
{
if(p[1]<=0) /* Point is in southern hemisphere, it */
{ /* doesn't matter if image passes equator! */
@@ -472,6 +501,7 @@ radecinimg(double *p, double *i, double *s, double *c)
}
}
}
+
return 0;
}
@@ -491,21 +521,25 @@ radecinimg(double *p, double *i, double *s, double *c)
there is an overlap (the survey image is completely within the crop
box). So we have to check both. */
int
-radecoverlap(struct onecropparams *crp)
+wcsmode_overlap(struct onecropparams *crp)
{
double *d, *fd;
double *i, *s, *c; /* for clear viewing. */
struct cropparams *p=crp->p;
+ size_t ndim=crp->p->imgs->ndim;
- /* First check if the four sides of the crop box are in the image.*/
- fd=(d=crp->corners)+8;
+ /* First check if the corners of the crop are in the image.*/
s=p->imgs[crp->in_ind].sized;
i=p->imgs[crp->in_ind].corners;
c=p->imgs[crp->in_ind].equatorcorr;
+ fd=(d=crp->corners) + 8;
do
{
- if( radecinimg(d, i, s, c) ) return 1;
- d+=2;
+ /* As long as one of the crop corners are in the image, we know there
+ is overlap and can return a true value. We don't need to check all
+ corners. */
+ if( point_in_dataset(d, i, s, c, ndim) ) return 1;
+ d+=ndim;
}
while(d<fd);
@@ -514,13 +548,14 @@ radecoverlap(struct onecropparams *crp)
s=crp->sized;
i=crp->corners;
c=crp->equatorcorr;
- fd=(d=p->imgs[crp->in_ind].corners)+8;
+ fd=(d=p->imgs[crp->in_ind].corners) + 8;
do
{
- if( radecinimg(d, i, s, c) ) return 1;
- d+=2;
+ if( point_in_dataset(d, i, s, c, ndim) ) return 1;
+ d+=ndim;
}
while(d<fd);
+ /* If control reaches here, there was no overlap. */
return 0;
}
diff --git a/bin/crop/wcsmode.h b/bin/crop/wcsmode.h
index 2866ddd..ead7fff 100644
--- a/bin/crop/wcsmode.h
+++ b/bin/crop/wcsmode.h
@@ -24,15 +24,15 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#define WCSMODE_H
void
-wcs_check_prepare(struct cropparams *p, struct inputimgs *img);
+wcsmode_check_prepare(struct cropparams *p, struct inputimgs *img);
void
-setcsides(struct onecropparams *crp);
+wcsmode_crop_corners(struct onecropparams *crp);
void
fillcrpipolygon(struct onecropparams *crp);
int
-radecoverlap(struct onecropparams *crp);
+wcsmode_overlap(struct onecropparams *crp);
#endif
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index d768b1d..ddcc89a 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -218,6 +218,7 @@ mkprof_build(void *inparam)
size_t i, id;
int lockresult;
+ double center[2];
long lpixel_o[2];
pthread_mutex_t *qlock=&p->qlock;
struct builtqueue *ibq, *fbq=NULL;
@@ -252,13 +253,15 @@ mkprof_build(void *inparam)
/* Get the overlapping pixels using the starting points (NOT
oversampled). */
- gal_box_border_from_center(p->x[id], p->y[id], mkp->width,
+ center[0]=p->x[id];
+ center[1]=p->y[id];
+ gal_box_border_from_center(center, 2, mkp->width,
ibq->fpixel_i, ibq->lpixel_i);
mkp->fpixel_i[0]=ibq->fpixel_i[0];
mkp->fpixel_i[1]=ibq->fpixel_i[1];
ibq->overlaps = gal_box_overlap(mkp->onaxes, ibq->fpixel_i,
ibq->lpixel_i, ibq->fpixel_o,
- lpixel_o);
+ lpixel_o, 2);
/* Build the profile if necessary, After this, the width is
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 2f7758a..bc36fdb 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -465,7 +465,7 @@ ui_matrix_make_align(struct warpparams *p, double *tmatrix)
/* Find the pixel scale along the two dimensions. Note that we will be
using the scale along the image X axis for both values. */
w=gal_wcs_warp_matrix(p->input->wcs);
- ps=gal_wcs_pixel_scale_deg(p->input->wcs);
+ ps=gal_wcs_pixel_scale(p->input->wcs);
/* Lets call the given WCS orientation `W', the rotation matrix we want
diff --git a/bin/warp/warp.c b/bin/warp/warp.c
index ff53451..4584ec4 100644
--- a/bin/warp/warp.c
+++ b/bin/warp/warp.c
@@ -455,7 +455,7 @@ correct_wcs_save_output(struct warpparams *p)
signs are usually different.*/
if( wcs->pc[1]<ABSOLUTEFLTERROR ) wcs->pc[1]=0.0f;
if( wcs->pc[2]<ABSOLUTEFLTERROR ) wcs->pc[2]=0.0f;
- pixelscale=gal_wcs_pixel_scale_deg(wcs);
+ pixelscale=gal_wcs_pixel_scale(wcs);
diff=fabs(wcs->pc[0])-fabs(wcs->pc[3]);
if( fabs(diff/pixelscale[0])<RELATIVEFLTERROR )
wcs->pc[3] = ( (wcs->pc[3] < 0.0f ? -1.0f : 1.0f) * fabs(wcs->pc[0]) );
diff --git a/configure.ac b/configure.ac
index df71c67..c738a63 100644
--- a/configure.ac
+++ b/configure.ac
@@ -50,7 +50,7 @@ AC_CONFIG_MACRO_DIRS([bootstrapped/m4])
# Library version, see the GNU Libtool manual ("Library interface versions"
# section for the exact definition of each) for
-GAL_CURRENT=1
+GAL_CURRENT=2
GAL_REVISION=0
GAL_AGE=0
GAL_LT_VERSION="${GAL_CURRENT}:${GAL_REVISION}:${GAL_AGE}"
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 077ebda..cccb440 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -2521,7 +2521,7 @@ you have downloaded
$ tar xf cfitsio_latest.tar.gz
$ cd cfitsio
$ ./configure --prefix=/usr/local --enable-sse2 --enable-reentrant
-$ make -j8 # Replace 8 with no. CPU threads.
+$ make
$ make utils
$ ./testprog > testprog.lis
$ diff testprog.lis testprog.out # Should have no output
@@ -2565,8 +2565,8 @@ downloaded
@url{ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib.tar.bz2,
@example
$ tar xf wcslib.tar.bz2
$ cd wcslib-X.X # Replace X.X with version number
-$ ./configure --without-pgplot LIBS="-pthread -lm"
-$ make -j8 # Replace 8 with no. CPU threads.
+$ ./configure --without-pgplot LIBS="-pthread -lm" --disable-fortran
+$ make
$ make check
$ sudo make install
@end example
@@ -7700,55 +7700,67 @@ image, see @ref{Crop section syntax}.
@node Crop modes, Crop section syntax, Crop, Crop
@subsection Crop modes
-In order to be comprehensive, intuitive, and easy to use, Crop has two ways
-to define crop region: 1) From its center and (square) side length, for
-example to generate postage stamps of a given catalog. 2) The vertices of
-the crop region, this can be useful for larger crops over many targets, for
-example to crop out a uniformly deep, or contiguous, region of a large
-survey. Irrespective of how the crop region is defined, both Image/pixel,
-and World coordinate system (WCS) coordinates are acceptable and all
-coordinates are read as floating point numbers (not integers, except for
-the @option{--section} option, see below). The coordinate standards used
-are the main modes of Crop. Here, the different ways to specify the crop
-region are discussed within each standard. For the full list options,
-please see @ref{Invoking astcrop}.
+In order to be comprehensive, intuitive, and easy to use, there are two
+ways to define the crop:
+
address@hidden
address@hidden
+From its center and side length. For example if you already know the
+coordinates of an object and want to inspect it in an image or to generate
+postage stamps of a catalog containing many such coordinates.
+
address@hidden
+The vertices of the crop region, this can be useful for larger crops over
+many targets, for example to crop out a uniformly deep, or contiguous,
+region of a large survey.
address@hidden enumerate
+
+Irrespective of how the crop region is defined, the coordinates to define
+the crop can be in Image (pixel) or World Coordinate System (WCS)
+standards. All coordinates are read as floating point numbers (not
+integers, except for the @option{--section} option, see below). By setting
+the @emph{mode} in Crop, you define the standard that the given coordinates
+must be interpretted. Here, the different ways to specify the crop region
+are discussed within each standard. For the full list options, please see
address@hidden astcrop}.
When the crop is defined by its center, the respective (integer) central
pixel position will be found internally according to the FITS standard. To
have this pixel positioned in the center of the cropped region, the final
-cropped region must have (in Image mode), or will have (in WCS mode) an add
-number of pixels. Furthermore, when the crop is defined as by its center,
-Crop allows you to only keep crops what don't have any blank pixels in the
-vicinity of their center (your primary target). This can be very convenient
-when your input catalog/coordinates originated from another survey/filter
-which is not fully covered by your input image, to learn more about this
-feature, please see the description of the @option{--checkcenter} option in
address@hidden astcrop}.
+cropped region will have an add number of pixels (even if you give an even
+number to @option{--width} in image mode).
+
+Furthermore, when the crop is defined as by its center, Crop allows you to
+only keep crops what don't have any blank pixels in the vicinity of their
+center (your primary target). This can be very convenient when your input
+catalog/coordinates originated from another survey/filter which is not
+fully covered by your input image, to learn more about this feature, please
+see the description of the @option{--checkcenter} option in @ref{Invoking
+astcrop}.
@table @asis
@item Image coordinates
-In image mode, Crop uses the pixel coordinates/positions (instead of world
-coordinates). In image mode, only one image may be input. The crop(s) can
-be defined in multiple ways as listed below.
+In image mode (@option{--mode=img}), Crop interprets the pixel coordinates
+and widths in units of the input data-elements (for example pixels in an
+image, not world coordinates). In image mode, only one image may be
+input. The output crop(s) can be defined in multiple ways as listed below.
@table @asis
@item Center of multiple crops (in a catalog)
-The center of (possibly multiple) crops are read from a text file. A
-catalog can also contain WCS coordinates, so @option{--mode=img} is
-necessary to avoid ambiguity when all columns are specified. The
address@hidden and @option{--ycol} columns in the catalog are interpreted
-as the center of a square crop with a width of @option{--iwidth} pixels (an
-odd number). The columns can contain any floating point value. The value to
address@hidden option is seen as a directory which will host (the
-possibly multiple) separate crop files, see @ref{Crop output} for more. For
-a tutorial using this feature, please see @ref{Hubble visually checks and
-classifies his catalog}.
+The center of (possibly multiple) crops are read from a text file. In this
+mode, the columns identified with the @option{--coordcol} option are
+interpreted as the center of a crop with a width of @option{--width} pixels
+along each dimension. The columns can contain any floating point value. The
+value to @option{--output} option is seen as a directory which will host
+(the possibly multiple) separate crop files, see @ref{Crop output} for
+more. For a tutorial using this feature, please see @ref{Hubble visually
+checks and classifies his catalog}.
@item Center of a single crop (on the command-line)
-The center of the crop is given on the command-line with the @option{--xc}
-and @option{--yc} options. The crop width is specified by the
address@hidden option. The given coordinates can be any floating point
-and the width has to be an odd integer, see the explanations above.
+The center of the crop is given on the command-line with the
address@hidden option. The crop width is specified by the
address@hidden option along each dimension. The given coordinates and
+width can be any floating point number.
@item Vertices of a single crop
In Image mode there are two options to define the vertices of a region to
@@ -7760,17 +7772,17 @@ section syntax} for a full description of this method.
The latter option (@option{--polygon}) is a higher-level method to define
any convex polygon (with any number of vertices) with floating point
values. Please see the description of this option in @ref{Invoking astcrop}
-for its syntax. It is available in both Image and WCS modes, hence, to read
-the vertices as pixel coordinates, @option{--mode=img} has to be called.
+for its syntax.
@end table
@item WCS coordinates
-Use the World Coordinate System (WCS) information to define the crop, not
-pixel coordinates. In this mode, the width (@option{--wwidth}) is read in
-units of arc-seconds and multiple images (tiles in a survey) can be
-input. When the cropped region (defined by center or vertices) overlaps
-with multiple of the input images/tiles, the overlapping regions will be
-taken from the respective input (they will be stitched in the crop).
+In WCS mode (@option{--mode=wcs}), the coordinates and widths are
+interpretted using the World Coordinate System (WCS, that must accompany
+the dataset), not pixel coordinates. In WCS mode, Crop accepts multiple
+datasets as input. When the cropped region (defined by its center or
+vertices) overlaps with multiple of the input images/tiles, the overlapping
+regions will be taken from the respective input (they will be stitched when
+necessary for each output crop).
In this mode, the input images do not necessarily have to be the same size,
they just need to have the same orientation and pixel resolution. Currently
@@ -7790,23 +7802,22 @@ listed below.
@table @asis
@item Center of multiple crops (in a catalog)
-Similar to catalog inputs in Image mode (above), but @option{--mode=wcs}
-has to be activated to avoid ambiguity. The central RA and Dec value for
-each crop will be read from the @option{--racol} and @option{--deccol}
-columns of the input catalog. The square cropped box will have an odd
-number of pixels.
+Similar to catalog inputs in Image mode (above), except that the values
+along each dimension are assumed to have the same units as the dataset's
+WCS information. For example, the central RA and Dec value for each crop
+will be read from the first and second calls to the @option{--coordcol}
+option. The width of the cropped box (in units of the WCS, or degrees in RA
+and Dec mode) must be specified with the @option{--width} option.
@item Center of a single crop (on the command-line)
-You can specify the center of only one crop box with the @option{--ra} and
address@hidden options. If it exists in the input images, it will be
-cropped similar to the catalog mode.
+You can specify the center of only one crop box with the @option{--center}
+option. If it exists in the input images, it will be cropped similar to the
+catalog mode, see above also for @code{--width}.
@item Vertices of a single crop
The @option{--polygon} option is a high-level method to define any convex
polygon (with any number of vertices). Please see the description of this
-option in @ref{Invoking astcrop} for its syntax. It is available in both
-Image and WCS modes, hence, to read the vertices as RA and Dec coordinates,
address@hidden has to be called/activated.
+option in @ref{Invoking astcrop} for its syntax.
@end table
@cartouche
@@ -7822,16 +7833,9 @@ input image with these coordinates, see @ref{Warp}.
@end table
As a summary, if you don't specify a catalog, you have to define the
-cropped region manually on the command-line. Other than the
address@hidden option, the other methods to define a single crop are
-unique to each mode, so the mode (@option{--mode}) will be ignored. When
-using a catalog to define the crop centers, if the columns to only one mode
-are given (for example only @option{--xcol} and @option{--ycol}, not the
-WCS columns) then the mode is irrelevant. However, when all image and WCS
-mode columns are given, @option{--mode} is mandatory (you can keep the
-columns in a configuration file and simply set the mode, see
address@hidden files} and @ref{Selecting table columns}). When using
address@hidden the mode is mandatory to interpret the values correctly.
+cropped region manually on the command-line. In any case the mode is
+mandatory for Crop to be able to interpret the values given as coordinates
+or widths.
@node Crop section syntax, Blank pixels, Crop modes, Crop
@@ -7942,32 +7946,32 @@ $ astcrop [OPTION...] [ASCIIcatalog] ASTRdata ...
One line examples:
@example
-## Crop all objects in catalog.txt from image.fits:
-$ astcrop catalog.txt image.fits
+## Crop all objects in cat.txt from image.fits:
+$ astcrop --catalog=cat.txt image.fits
## Crop all options in catalog (with RA,DEC) from all the files
## ending in `_drz.fits' in `/mnt/data/COSMOS/':
-$ astcrop --mode=wcs catalog.txt /mnt/data/COSMOS/*_drz.fits
+$ astcrop --mode=wcs --catalog=cat.txt /mnt/data/COSMOS/*_drz.fits
-## Crop-out the outer 10 border pixels of input image:
+## Crop the outer 10 border pixels of the input image:
$ astcrop --section=10:*-10,10:*-10 --hdu=2 image.fits
## Crop region around RA and Dec of (189.16704, 62.218203):
-$ astcrop --ra=189.16704 --dec=62.218203 goodsnorth.fits
+$ astcrop --mode=wcs --center=189.16704,62.218203 goodsnorth.fits
## Crop region around pixel coordinate (568.342, 2091.719):
-$ astcrop --xc=568.342 --yc=2091.719 --iwidth=201 image.fits
+$ astcrop --mode=img --center=568.342,2091.719 --width=201 image.fits
@end example
@noindent
Crop has one mandatory argument which is the input image name(s), shown
above with @file{ASTRdata ...}. You can use shell expansions, for example
@command{*} for this if you have lots of images in WCS mode. If the crop
-box centers are in a catalog, you also have to provide the catalog name as
-an argument. Alternatively, you have to provide the crop box parameters
-with command-line options. See @ref{Crop output} for how the output file
-name(s) can be specified. For the full list of general options to all
-Gnuastro programs (including Crop), please see @ref{Common options}.
+box centers are in a catalog, you can use the @option{--catalog} option. In
+other cases, you have to provide the single cropped output parameters must
+be given with command-line options. See @ref{Crop output} for how the
+output file name(s) can be specified. For the full list of general options
+to all Gnuastro programs (including Crop), please see @ref{Common options}.
Floating point numbers can be used to specify the crop region (except the
@option{--section} option, see @ref{Crop section syntax}). In such cases,
@@ -7975,11 +7979,10 @@ the floating point values will be used to find the
desired integer pixel
indices based on the FITS standard. Hence, Crop ultimately doesn't do any
sub-pixel cropping (in other words, it doesn't change pixel values). If you
need such crops, you can use @ref{Warp} to first warp the image to the a
-new pixel grid where your initial floating points can be seen as integers,
-then crop from that with Crop. For example, let's say you want a crop from
-pixels 12.982 to 80.982 along the first dimension. You should first
-translate the image by @mymath{-0.482} (note that the edge of a pixel is at
-integer multiples of @mymath{0.5}). So you should run Warp with
+new pixel grid, then crop from that. For example, let's assume you want a
+crop from pixels 12.982 to 80.982 along the first dimension. You should
+first translate the image by @mymath{-0.482} (note that the edge of a pixel
+is at integer multiples of @mymath{0.5}). So you should run Warp with
@option{--translate=-0.482,0} and then crop the warped image with
@option{--section=13:81}.
@@ -7992,12 +7995,12 @@ for more.
@cindex Asynchronous thread allocation
When in catalog mode, Crop will run in parallel unless you set
address@hidden, see @ref{Multi-threaded operations}. Note that when
-multiple threads are being used, in verbose mode, the outputs will not be
-in order. This is because the threads are asynchronous and thus not started
-in order. This has no effect on the output, see @ref{Hubble visually checks
-and classifies his catalog} for a tutorial on effectively using this
-feature.
address@hidden, see @ref{Multi-threaded operations}. Note that
+when multiple outputs are created with threads, the outputs will not be
+created in the same order. This is because the threads are asynchronous and
+thus not started in order. This has no effect on each output, see
address@hidden visually checks and classifies his catalog} for a tutorial on
+effectively using this feature.
@menu
* Crop options:: A list of all the options with explanation.
@@ -8021,14 +8024,6 @@ example, when using @option{--section}, the right
ascension of the bottom
left and top left corners will not be equal. If you only want regions
within a given right ascension, use @option{--polygon} in WCS mode.
address@hidden
address@hidden
address@hidden:} The coordinates are in the FITS format. So the first
-axis is the horizontal axis when viewed in SAO ds9 and the second axis
-is the vertical. Also in the FITS standard, counting begins from 1
-(one) not 0 (zero).
address@hidden cartouche
-
@noindent
Input image parameters:
@table @option
@@ -8064,42 +8059,33 @@ coordinate system on the input images. See
@option{--hstartwcs}
Crop box parameters:
@table @option
address@hidden -x FLT
address@hidden --xc=FLT
-X axis (first image axis, horizontal when viewed in SAO ds9) position of
-crop box center. Along with @option{--yc}, this only produces one crop in
-each run. The width of the square crop (in pixels) is the value to
address@hidden
-
address@hidden -y FLT
address@hidden --yc=FLT
-Y axis (second image axis, vertical when viewed in SAO ds9) position of
-crop box center. See @option{--xc}.
-
address@hidden -r FLT
address@hidden --ra=FLT
-Right Ascension of crop box center. Along with @option{--dec}, this only
-produces one crop in each run. The width of the square crop (in arcseconds)
-is the value to @option{--wwidth}.
-
address@hidden -d FLT
address@hidden --dec=FLT
-Declination of crop box center. see @option{--ra}.
-
address@hidden -a INT
address@hidden --iwidth=INT
-Width of the cropped region in units of pixels when the crop is defined by
-its center in Image mode (see @ref{Crop modes}). In order for the chosen
-central pixel to be in the center of the cropped image, this option only
-accepts an odd number (an error will be given if its even). If you want an
-even sided crop, you can run Crop afterwards with
address@hidden -c FLT[,FLT[,...]]
address@hidden --center=FLT[,FLT[,...]]
+The central position of the crop in the input image. The positions along
+each dimension must be separated by a comma (@key{,}) and fractions are
+also acceptable. The number of values given to this option must be the same
+as the dimensions of the input dataset. The width of the crop should be set
+with @code{--width}. The units of the coordinates are read based on the
+value to the @option{--mode} option, see below.
+
address@hidden -w FLT[,FLT[,...]]
address@hidden --width=FLT[,FLT[,...]]
+Width of the cropped region about its center. @option{--width} may take
+either a single value (to be used for all dimensions) or multiple values (a
+specific value for each dimension). If in WCS mode, value(s) given to this
+option will be read in the same units as the dataset's WCS information
+along this dimension. The final output will have an odd number of pixels to
+allow easy identification of the pixel which keeps your requested
+coordinate (from @option{--center} or @option{--catalog}).
+
+The @code{--width} option also accepts fractions. For example if you want
+the width of your crop to be 3 by 5 arcseconds along RA and Dec
+respectively, you can call it with: @option{--width=3/3600,5/3600}.
+
+If you want an even sided crop, you can run Crop afterwards with
@option{--section=":*-1,:*-1"} or @option{--section=2:,2:} (depending on
which side you don't need), see @ref{Crop section syntax}.
address@hidden -w FLT
address@hidden --wwidth=FLT
-The width of the crop box in WCS mode in units of arc-seconds.
-
@item -l STR
@itemx --polygon=STR
String of crop polygon vertices. Note that currently only convex polygons
@@ -8195,31 +8181,12 @@ Section of the input image which you want to be
cropped. See @ref{Crop
section syntax} for a complete explanation on the syntax required for this
input.
address@hidden -i STR/INT
address@hidden --xcol=STR/INT
-Column selection for the first FITS axis position of the box center in a
-catalog. In SAO ds9, the first FITS axis is the horizontal axis. The value
-can be either the column number (starting from 1), or a match/search in the
-table meta-data, see @ref{Selecting table columns}.
-
address@hidden -j STR/INT
address@hidden --ycol=STR/INT
-Column selection for the second FITS axis position of the box center. In
-SAO ds9, the second FITS axis is the vertical axis. The value can be either
-the column number (starting from 1), or a match/search in the table
-meta-data, see @ref{Selecting table columns}.
-
address@hidden -f STR/INT
address@hidden --racol=STR/INT
-Column selection for Right Ascension (RA) in the input catalog. The value
-can be either the column number (starting from 1), or a match/search in the
-table meta-data, see @ref{Selecting table columns}.
-
address@hidden -g STR/INT
address@hidden --deccol=STR/INT
-Column selection of declination in the input catalog. The value can be
-either the column number (starting from 1), or a match/search in the table
-meta-data, see @ref{Selecting table columns}.
address@hidden -x STR/INT
address@hidden --coordcol=STR/INT
+The column in a catalog to read as a coordinate. This option must be called
+multiple times depending on the number of dimensions in the input
+dataset. The value can be either the column number (starting from 1), or a
+match/search in the table meta-data, see @ref{Selecting table columns}.
@item -n STR/INT
@item --namecol=STR/INT
@@ -8308,18 +8275,21 @@ image or WCS. The value must either be @option{img} or
@option{wcs}, see
@node Crop output, , Crop options, Invoking astcrop
@subsubsection Crop output
-The value given to @option{--output} option will depend on how many crops
-were created, see @ref{Crop modes}:
+The string given to @option{--output} option will be interpretted depending
+on how many crops were requested, see @ref{Crop modes}:
@itemize
@item
When a catalog is given, the value of the @option{--output} (see
@ref{Common options}) will be read as the directory to store the output
cropped images. Hence if it doesn't already exist, Crop will abort with an
-error of a ``No such file or directory'' error. The crop file names will
-consist of two parts: a variable part (the row number of each target
-starting from 1) along with a fixed string which you can set with the
address@hidden option.
+error of a ``No such file or directory'' error.
+
+The crop file names will consist of two parts: a variable part (the row
+number of each target starting from 1) along with a fixed string which you
+can set with the @option{--suffix} option. Optionally, you may also use the
address@hidden option to define a column in the input catalog to use as
+the file name instead of numbers.
@item
When only one crop is desired, the value to @option{--output} will be read
@@ -8341,9 +8311,11 @@ range from that input image in the same syntax as
@ref{Crop section
syntax}. So this string can be directly given to the @option{--section}
option later.
-Once done, a log file will be created in the current directory named
address@hidden This file will have three columns and the same
-number of rows as the number of cropped images.
+Once done, a log file can be created in the current directory with the
address@hidden option. This file will have three columns and the same number
+of rows as the number of cropped images. There are also comments on the top
+of the log file explaining basic information about the run and descriptions
+for the columns. A short description of the columns is also given below:
@enumerate
@item
@@ -8358,14 +8330,6 @@ given a value of 0 (see @ref{Invoking astcrop}), the
center will not be
checked and this column will be given a value of @code{-1}.
@end enumerate
-There are also comments on the top of the log file explaining basic
-information about the run and descriptions for the columns. If the log file
-cannot be created (for example you don't have write permission in the
-directory you are running Crop in) or you have specifically asked for no
-log file (with the @option{--nolog} option), then a log file will not be
-created (unless @option{--individual} is called).
-
-
@@ -19669,7 +19633,7 @@ which is better considering floating point errors:
@dispmath{{\sin^2(d)\over 2}=\sin^2\left( {d_1-d_2\over 2}
\right)+\cos(d_1)\cos(d_2)\sin^2\left( {r_1-r_2\over 2} \right)}
@end deftypefun
address@hidden {double *} gal_wcs_pixel_scale_deg (struct wcsprm @code{*wcs})
address@hidden {double *} gal_wcs_pixel_scale (struct wcsprm @code{*wcs})
Return the pixel scale for each dimension of @code{wcs} in degrees. The
output is an array of double precision floating point type with one element
for each dimension.
diff --git a/lib/blank.c b/lib/blank.c
index 487f77f..44dee32 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -195,7 +195,7 @@ gal_blank_present(gal_data_t *input, int updateflag)
error(EXIT_FAILURE, 0, "%s: tile mode is currently not supported for "
"strings", __func__);
strf = (str=input->array) + input->size;
- do if(!strcmp(*str++,GAL_BLANK_STRING)) return 1; while(str<strf);
+ do if(!strcmp(*str,GAL_BLANK_STRING)) return 1; while(++str<strf);
break;
/* Complex types */
diff --git a/lib/box.c b/lib/box.c
index 365d027..d7d3386 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -24,6 +24,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <math.h>
#include <stdio.h>
+#include <errno.h>
+#include <error.h>
#include <stdlib.h>
#include <gnuastro/box.h>
@@ -86,29 +88,40 @@ gal_box_ellipse_in_box(double a, double b, double
theta_rad, long *width)
/* We have the central pixel and box size of the crop box, find the
- starting and ending pixels: */
+ starting (first) and ending (last) pixels: */
void
-gal_box_border_from_center(double xc, double yc, long *width,
+gal_box_border_from_center(double *center, size_t ndim, long *width,
long *fpixel, long *lpixel)
{
- long lxc, lyc;
+ size_t i;
+ long tmp;
double intpart;
- /* Round the double values in a the long values: */
- lxc=xc;
- if (fabs(modf(xc, &intpart))>=0.5)
- ++lxc;
- lyc=yc;
- if (fabs(modf(yc, &intpart))>=0.5)
- ++lyc;
-
- /* Set the initial values for the actual image: */
- fpixel[0]=lxc-width[0]/2; fpixel[1]=lyc-width[1]/2;
- lpixel[0]=lxc+width[0]/2; lpixel[1]=lyc+width[1]/2;
- /*
- printf("\n\nCenter is on: %ld, %ld\n", lxc, lyc);
- printf("Starting and ending pixels: (%ld, %ld) -- (%ld, %ld)\n\n\n",
- fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
+ /* Go over the dimensions. */
+ for(i=0;i<ndim;++i)
+ {
+ /* When the decimal fraction of the center (in floating point) is
+ larger than 0.5, then it is must be incremented. */
+ tmp = center[i] + ( fabs(modf(center[i], &intpart))>=0.5 ? 1 : 0 );
+
+ /* Set the first and last pixels in this dimension. */
+ fpixel[i]=tmp-width[i]/2;
+ lpixel[i]=tmp+width[i]/2;
+ }
+
+ /* For a check:
+ if(ndim==2)
+ {
+ printf("center: %g, %g\n", center[0], center[1]);
+ printf("fpixel: %ld, %ld\n", fpixel[0], fpixel[1]);
+ printf("lpixel: %ld, %ld\n", lpixel[0], lpixel[1]);
+ }
+ else if(ndim==3)
+ {
+ printf("center: %g, %g, %g\n", center[0], center[1], center[2]);
+ printf("fpixel: %ld, %ld, %ld\n", fpixel[0], fpixel[1], fpixel[2]);
+ printf("lpixel: %ld, %ld, %ld\n", lpixel[0], lpixel[1], lpixel[2]);
+ }
*/
}
@@ -116,6 +129,7 @@ gal_box_border_from_center(double xc, double yc, long
*width,
+
/* Problem to solve: We have set the first and last pixels in an input
image (fpixel_i[2] and lpixel_i[2]). But those first and last
pixels don't necessarily have to lie within the image, they can be
@@ -167,7 +181,7 @@ gal_box_border_from_center(double xc, double yc, long
*width,
-----------------------------
- So in short the arrays we are dealing with here are:
+ So, in short the arrays we are dealing with here are:
fpixel_i: Coordinates of the first pixel in input image.
lpixel_i: Coordinates of the last pixel in input image.
@@ -188,80 +202,85 @@ gal_box_border_from_center(double xc, double yc, long
*width,
*/
int
gal_box_overlap(long *naxes, long *fpixel_i, long *lpixel_i,
- long *fpixel_o, long *lpixel_o)
+ long *fpixel_o, long *lpixel_o, size_t ndim)
{
- long width[2];
+ size_t i;
+ long width;
/* In case you want to see how things are going:
- printf("\n\nImage size: [%ld, %ld]\n", naxes[0], naxes[1]);
- printf("fpixel_i -- lpixel_i: (%ld, %ld) -- (%ld, %ld)\n\n", fpixel_i[0],
- fpixel_i[1], lpixel_i[0], lpixel_i[1]);
+ printf("\n\nImage size: [");
+ for(i=0;i<ndim;++i) {printf("%ld, ", naxes[i]);} printf("\b\b]\n");
+ printf("fpixel_i -- lpixel_i: (");
+ for(i=0;i<ndim;++i) {printf("%ld, ", fpixel_i[i]);} printf("\b\b) -- (");
+ for(i=0;i<ndim;++i) {printf("%ld, ", lpixel_i[i]);} printf("\b\b)\n");
*/
- width[0]=lpixel_i[0]-fpixel_i[0]+1;
- width[1]=lpixel_i[1]-fpixel_i[1]+1;
-
- /* Set the initial values for the cropped array: */
- fpixel_o[0]=1; fpixel_o[1]=1;
- lpixel_o[0]=width[0]; lpixel_o[1]=width[1];
-
-
- /* Check the four corners to see if they should be adjusted: To
- understand the first, look at this, suppose | separates pixels
- and the * shows where the image actually begins.
-
- |-2|-1| 0* 1| 2| 3| 4| (survey image)
- *1 | 2| 3| 4| 5| 6| 7| (crop image)
-
- So when fpixel_i is negative, e.g., fpixel_i=-2, then the index
- of the pixel in the cropped image we want to begin with, that
- corresponds to 1 in the survey image is: fpixel_o= 4 = 2 +
- -1*fpixel_i.*/
- if(fpixel_i[0]<1)
- {
- fpixel_o[0]=-1*fpixel_i[0]+2;
- fpixel_i[0]=1;
- }
- if(fpixel_i[1]<1)
+ /* Do the calculations for each dimension: */
+ for(i=0;i<ndim;++i)
{
- fpixel_o[1]=-1*fpixel_i[1]+2;
- fpixel_i[1]=1;
+ /* Set the width and initialize the overlap values. */
+ fpixel_o[i] = 1;
+ lpixel_o[i] = width = lpixel_i[i] - fpixel_i[i] + 1;
+
+
+ /* Check the four corners to see if they should be adjusted: To
+ understand the first, look at this, suppose | separates pixels and
+ the * shows where the image actually begins.
+
+ |-2|-1| 0* 1| 2| 3| 4| ( input image )
+ *1 | 2| 3| 4| 5| 6| 7| ( crop image )
+
+ So when fpixel_i is negative, e.g., fpixel_i=-2, then the index of
+ the pixel in the cropped image we want to begin with, that
+ corresponds to 1 in the survey image is:
+
+ fpixel_o = 4 = 2 + -1*fpixel_i
+ */
+ if(fpixel_i[i]<1)
+ {
+ fpixel_o[i] = -1*fpixel_i[i]+2;
+ fpixel_i[i] = 1;
+ }
+
+
+ /*The same principle applies to the end of an image. Take "s" is the
+ maximum size along a specific axis in the survey image and "c" is
+ the size along the same axis on the cropped image. Assume the the
+ cropped region's last pixel in that axis will be 2 pixels larger
+ than s:
+
+ |s-1| s* s+1| s+2| (survey image)
+ |c-3| c-2| c-1| c* (crop image)
+
+ So you see that if the outer pixel is n pixels away then in the
+ cropped image we should only fill upto c-n.*/
+ if(lpixel_i[i]>naxes[i])
+ {
+ lpixel_o[i] = width - (lpixel_i[i]-naxes[i]);
+ lpixel_i[i] = naxes[i];
+ }
}
+ /* In case you want to see the results.
+ printf("\nAfter correction:\n");
+ printf("Input image: (");
+ for(i=0;i<ndim;++i) {printf("%ld, ", fpixel_i[i]);} printf("\b\b) -- (");
+ for(i=0;i<ndim;++i) {printf("%ld, ", lpixel_i[i]);} printf("\b\b)\n");
+ printf("output image: (");
+ for(i=0;i<ndim;++i) {printf("%ld, ", fpixel_o[i]);} printf("\b\b) -- (");
+ for(i=0;i<ndim;++i) {printf("%ld, ", lpixel_o[i]);} printf("\b\b)\n");
+ */
- /*The same principle applies to the end of an image. Take "s"
- is the maximum size along a specific axis in the survey image
- and "c" is the size along the same axis on the cropped image.
- Assume the the cropped region's last pixel in that axis will
- be 2 pixels larger than s:
-
- |s-1| s* s+1| s+2| (survey image)
- |c-3| c-2| c-1| c* (crop image)
- So you see that if the outer pixel is n pixels away then in
- the cropped image we should only fill upto c-n.*/
- if (lpixel_i[0]>naxes[0])
- {
- lpixel_o[0]=width[0]-(lpixel_i[0]-naxes[0]);
- lpixel_i[0]=naxes[0];
- }
- if (lpixel_i[1]>naxes[1])
- {
- lpixel_o[1]=width[1]-(lpixel_i[1]-naxes[1]);
- lpixel_i[1]=naxes[1];
- }
+ /* If the first image pixel in every dimension is larger than the input
+ image's width or the last image pixel is before the input image, then
+ there is no overlap and we must return 0. */
+ for(i=0;i<ndim;++i)
+ if( fpixel_i[i]>naxes[i] || lpixel_i[i]<1 )
+ return 0;
- /* In case you wish to see the results.
- printf("\nAfter correction:\n");
- printf("Input image: (%ld, %ld) -- (%ld, %ld)\n", fpixel_i[0],
- fpixel_i[1], lpixel_i[0], lpixel_i[1]);
- printf("output image:(%ld, %ld) -- (%ld, %ld)\n\n\n", fpixel_o[0],
- fpixel_o[1], lpixel_o[0], lpixel_o[1]);
- */
- if(fpixel_i[0]>naxes[0] || fpixel_i[1]>naxes[1]
- || lpixel_i[0]<1 || lpixel_i[1]<1)
- return 0;
- else
- return 1;
+ /* The first and last image pixels are within the image, so there is
+ overlap. */
+ return 1;
}
diff --git a/lib/gnuastro-internal/options.h b/lib/gnuastro-internal/options.h
index caddc38..9a9f974 100644
--- a/lib/gnuastro-internal/options.h
+++ b/lib/gnuastro-internal/options.h
@@ -48,6 +48,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
only applies to them. */
#define GAL_OPTIONS_MAX_VALUE_LEN 10
+/* Statically allocated space for printing option values. */
+#define GAL_OPTIONS_STATIC_MEM_FOR_VALUES 2000
+
diff --git a/lib/gnuastro/box.h b/lib/gnuastro/box.h
index 5ee522a..a1aa764 100644
--- a/lib/gnuastro/box.h
+++ b/lib/gnuastro/box.h
@@ -55,12 +55,12 @@ void
gal_box_ellipse_in_box(double a, double b, double theta_rad, long *width);
void
-gal_box_border_from_center(double xc, double yc, long *width,
+gal_box_border_from_center(double *center, size_t ndim, long *width,
long *fpixel, long *lpixel);
int
gal_box_overlap(long *naxes, long *fpixel_i, long *lpixel_i,
- long *fpixel_o, long *lpixel_o);
+ long *fpixel_o, long *lpixel_o, size_t ndim);
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 226d534..314891e 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -65,8 +65,8 @@ gal_dimension_num_neighbors(size_t ndim);
/************************************************************************/
/******************** Coordinates **********************/
/************************************************************************/
-#define GAL_DIMENSION_FLT_TO_INT(FLT) ( (FLT)-(int)(FLT) > 0.5f \
- ? (int)(FLT)+1 : (int)(FLT) )
+#define GAL_DIMENSION_FLT_TO_INT(FLT) ( (FLT)-(long)(FLT) > 0.5f \
+ ? (long)(FLT)+1 : (long)(FLT) )
void
gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 8a2a70e..78c6dec 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -86,7 +86,7 @@ double
gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
double *
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs);
+gal_wcs_pixel_scale(struct wcsprm *wcs);
double
gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
diff --git a/lib/options.c b/lib/options.c
index bc57997..3c25003 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -514,7 +514,6 @@ gal_options_read_tableformat(struct argp_option *option,
char *arg,
element of the option. The last element of the array is `-1' to allow
finding the number of elements within it later (similar to a string
which terminates with a '\0' element). */
-#define PARSE_SIZES_STATICSTR_LEN 2000
void *
gal_options_parse_sizes_reverse(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk)
@@ -523,7 +522,7 @@ gal_options_parse_sizes_reverse(struct argp_option *option,
char *arg,
double *v;
gal_data_t *values;
size_t nc, num, *array;
- char *str, sstr[PARSE_SIZES_STATICSTR_LEN];
+ char *str, sstr[GAL_OPTIONS_STATIC_MEM_FOR_VALUES];
/* We want to print the stored values. */
if(lineno==-1)
@@ -537,12 +536,12 @@ gal_options_parse_sizes_reverse(struct argp_option
*option, char *arg,
nc=0;
for(i=num-1;i>=0;--i)
{
- if( nc > PARSE_SIZES_STATICSTR_LEN-100 )
+ if( nc > GAL_OPTIONS_STATIC_MEM_FOR_VALUES-100 )
error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so we "
"can address the problem. The number of necessary "
"characters in the statically allocated string has become "
"too close to %d", __func__, PACKAGE_BUGREPORT,
- PARSE_SIZES_STATICSTR_LEN);
+ GAL_OPTIONS_STATIC_MEM_FOR_VALUES);
nc += sprintf(sstr+nc, "%zu,", array[i]);
}
sstr[nc-1]='\0';
@@ -567,14 +566,17 @@ gal_options_parse_sizes_reverse(struct argp_option
*option, char *arg,
for(i=0;i<values->size;++i)
{
if(v[i]<0)
- error(EXIT_FAILURE, 0, "the given value in `%s' (%g) is not 0 "
- "or positive. The values to the `--%s' option must be "
- "positive", arg, v[i], option->name);
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the given "
+ "value in `%s' (%g) is not 0 or positive. The "
+ "values to the `--%s' option must be positive",
+ arg, v[i], option->name);
+
if(ceil(v[i]) != v[i])
- error(EXIT_FAILURE, 0, "the given value in `%s' (%g) is not an "
- "integer. The values to the `--%s' option must be "
- "integers", arg, v[i], option->name);
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the given "
+ "value in `%s' (%g) is not an integer. The "
+ "values to the `--%s' option must be integers",
+ arg, v[i], option->name);
}
/* Write the values into an allocated size_t array and finish it with
diff --git a/lib/wcs.c b/lib/wcs.c
index f99eea0..289dd52 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -363,7 +363,7 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
if(wcs->pc)
{
/* Get the pixel scale. */
- ps=gal_wcs_pixel_scale_deg(wcs);
+ ps=gal_wcs_pixel_scale(wcs);
/* The PC matrix and the CDELT elements might both contain scale
elements (during processing the scalings might be added only to PC
@@ -431,9 +431,9 @@ gal_wcs_angular_distance_deg(double r1, double d1, double
r2, double d2)
-/* Return the pixel scale of the image in both dimentions in degrees. */
+/* Return the pixel scale of the dataset in units of the WCS. */
double *
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs)
+gal_wcs_pixel_scale(struct wcsprm *wcs)
{
gsl_vector S;
gsl_matrix A, V;
@@ -490,7 +490,7 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
"The input WCS has %d dimensions", __func__, wcs->naxis);
/* Get the pixel scales along each axis in degrees, then multiply. */
- pixscale=gal_wcs_pixel_scale_deg(wcs);
+ pixscale=gal_wcs_pixel_scale(wcs);
/* Clean up and return the result. */
out = pixscale[0] * pixscale[1] * 3600.0f * 3600.0f;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 73190ca..3948ce7 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -87,18 +87,18 @@ if COND_COSMICCAL
cosmiccal/simpletest.sh: prepconf.sh.log
endif
if COND_CROP
- MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/xcyc.sh
\
- crop/xcycnoblank.sh crop/section.sh crop/radec.sh crop/imgpolygon.sh \
- crop/imgoutpolygon.sh crop/wcspolygon.sh
+ MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/imgcenter.sh \
+ crop/imgcenternoblank.sh crop/section.sh crop/wcscenter.sh \
+ crop/imgpolygon.sh crop/imgoutpolygon.sh crop/wcspolygon.sh
crop/imgcat.sh: mkprof/mosaic1.sh.log
crop/wcscat.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
- crop/xcyc.sh: mkprof/mosaic1.sh.log
- crop/xcycnoblank.sh: mkprof/mosaic1.sh.log
+ crop/imgcenter.sh: mkprof/mosaic1.sh.log
+ crop/imgcenternoblank.sh: mkprof/mosaic1.sh.log
crop/section.sh: mkprof/mosaic1.sh.log
- crop/radec.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
- mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
+ crop/wcscenter.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
+ mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
crop/imgpolygon.sh: mkprof/mosaic1.sh.log
crop/imgoutpolygon.sh: mkprof/mosaic1.sh.log
crop/wcspolygon.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
diff --git a/tests/crop/imgcat.sh b/tests/crop/imgcat.sh
index b782fcd..824db0d 100755
--- a/tests/crop/imgcat.sh
+++ b/tests/crop/imgcat.sh
@@ -54,5 +54,6 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# enable multithreaded access to files, the tests pass. It is the
# users choice to enable this feature.
cat=$topsrc/tests/$prog/cat.txt
-$execname $img --catalog=$cat --mode=img --suffix=_imgcat.fits --numthreads=1 \
- --zeroisnotblank --xcol=X_CENTER --ycol=Y_CENTER --namecol=NAME
+$execname $img --catalog=$cat --suffix=_imgcat.fits --numthreads=1 \
+ --zeroisnotblank --coordcol=X_CENTER --coordcol=Y_CENTER \
+ --namecol=NAME --mode=img --width=201
diff --git a/tests/crop/xcyc.sh b/tests/crop/imgcenter.sh
similarity index 89%
rename from tests/crop/xcyc.sh
rename to tests/crop/imgcenter.sh
index f472ee0..f56d149 100755
--- a/tests/crop/xcyc.sh
+++ b/tests/crop/imgcenter.sh
@@ -1,4 +1,4 @@
-# Crop a box based on --xc and --yc.
+# Make a single crop with center defined in Image mode.
#
# See the Tests subsection of the manual for a complete explanation
# (in the Installing gnuastro section).
@@ -53,4 +53,5 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# The number of threads is one so if CFITSIO does is not configured to
# enable multithreaded access to files, the tests pass. It is the
# users choice to enable this feature.
-$execname $img --xc=251 --yc=251 --output=crop_xcyc.fits --numthreads=1
+$execname $img --center=251,251 --output=crop_imgcenter.fits \
+ --numthreads=1 --mode=img --width=201
diff --git a/tests/crop/xcycnoblank.sh b/tests/crop/imgcenternoblank.sh
similarity index 86%
rename from tests/crop/xcycnoblank.sh
rename to tests/crop/imgcenternoblank.sh
index 394ab4a..e23826b 100755
--- a/tests/crop/xcycnoblank.sh
+++ b/tests/crop/imgcenternoblank.sh
@@ -1,4 +1,5 @@
-# Crop a box based on --xc and --yc, with no blank pixels.
+# Crop a region defined by its center in image coordinates, with no blank
+# pixels.
#
# See the Tests subsection of the manual for a complete explanation
# (in the Installing gnuastro section).
@@ -53,5 +54,5 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# The number of threads is one so if CFITSIO does is not configured to
# enable multithreaded access to files, the tests pass. It is the
# users choice to enable this feature.
-$execname $img --xc=500 --yc=500 --noblank --output=crop_xcycnb.fits \
- --numthreads=1
+$execname $img --center=500,500 --noblank --numthreads=1 \
+ --output=crop_imgcenternoblank.fits --mode=img --width=201
diff --git a/tests/crop/section.sh b/tests/crop/section.sh
index 5df61bf..efcc623 100755
--- a/tests/crop/section.sh
+++ b/tests/crop/section.sh
@@ -53,4 +53,4 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# enable multithreaded access to files, the tests pass. It is the
# users choice to enable this feature.
$execname $img --section=-10:*+10,:250 --output=crop_section.fits \
- --numthreads=1
+ --numthreads=1 --mode=img
diff --git a/tests/crop/wcscat.sh b/tests/crop/wcscat.sh
index 011194b..48a1923 100755
--- a/tests/crop/wcscat.sh
+++ b/tests/crop/wcscat.sh
@@ -55,5 +55,6 @@ done
# enable multithreaded access to files, the tests pass. It is the
# users choice to enable this feature.
cat=$topsrc/tests/$prog/cat.txt
-$execname $img --catalog=$cat --mode=wcs --suffix=_wcscat.fits \
- --zeroisnotblank --racol=4 --deccol=DEC_CENTER --numthreads=1
+$execname $img --catalog=$cat --suffix=_wcscat.fits \
+ --zeroisnotblank --coordcol=4 --coordcol=DEC_CENTER \
+ --numthreads=1 --mode=wcs --width=3/3600
diff --git a/tests/crop/radec.sh b/tests/crop/wcscenter.sh
similarity index 88%
rename from tests/crop/radec.sh
rename to tests/crop/wcscenter.sh
index 001f961..5d62411 100755
--- a/tests/crop/radec.sh
+++ b/tests/crop/wcscenter.sh
@@ -1,4 +1,4 @@
-# Crop a box based on the --ra and --dec options.
+# Crop a box based on center in WCS coordinates.
#
# See the Tests subsection of the manual for a complete explanation
# (in the Installing gnuastro section).
@@ -49,5 +49,5 @@ done
# Actual test script
# ==================
-$execname $img --ra=0.99917157 --dec=1.0008283 --wwidth=0.3 \
- --output=crop_radec.fits
+$execname $img --center=0.99917157,1.0008283 --output=crop_wcscenter.fits \
+ --mode=wcs --width=0.3/3600
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index 36f6f26..578155c 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -70,7 +70,7 @@
# space characters in them, quote the full value
numjobs=8
builddir=build
-outdir=
+outdir=~/desktop
# Set the utility name, along with its arguments and options. NOTE, for
@@ -80,9 +80,9 @@ outdir=
# script, and once for the utility. In such cases it might be easier to
# just add the argument/option to the final script that runs the utility
# rather than these variables.
-utilname=
-arguments=
-options=
+utilname=crop
+arguments=/home/mohammad/documents/personal/research/software/local/c/gnuastro/build/tests/mkprofcat1.fits
+options="--section=-10:*+10,:250 --output=crop_section.fits --numthreads=1
--mode=img -h0"
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 58ec3af: Unified options used for single crop from center,
Mohammad Akhlaghi <=