/*********************************************************************
Crop - Crop a given size from one or multiple images.
Crop is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi
Contributing author(s):
Copyright (C) 2015-2018, Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
Gnuastro is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with Gnuastro. If not, see .
**********************************************************************/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "main.h"
#include "onecrop.h"
/*******************************************************************/
/**************** Check for ui.c *********************/
/*******************************************************************/
/* This function is called from ui.c. Its job is to check the WCS
values of this */
void
wcsmode_check_prepare(struct cropparams *p, struct inputimgs *img)
{
double *pixscale;
struct wcsprm *wcs=img->wcs;
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];
/* Check if the image is aligned with the WCS coordinates. Note that
because of small floating point errors, some programs might still keep
very small values in the off-diagonal matrix elements. */
if( fabs(wcs->pc[1]/wcs->pc[3])>1e-6 || fabs(wcs->pc[2]/wcs->pc[3])>1e-6 )
error(EXIT_FAILURE, 0, "%s: HDU %s: is not aligned to the "
"celestial coordinates. The first FITS axis should be "
"along the Right Ascension and the second FITS axis "
"should be along the declination.\n\n"
"Gnuastro's Warp program can align it with the following "
"command:\n\n"
" $ astwarp %s --hdu=%s --align\n",
img->name, p->cp.hdu, img->name, p->cp.hdu);
if(wcs->pc[0]>0)
error(EXIT_FAILURE, 0, "%s: HDU %s: An increase in the first "
"FITS axis pixel coordinates should be a decrease in the "
"RA. You have to flip the image along the second axis "
"before running Crop", img->name, p->cp.hdu);
if(wcs->pc[3]<0)
error(EXIT_FAILURE, 0, "%s: HDU %s: An increase in the second "
"FITS axis pixel coordinates should translate to an "
"increase in the declination. You have to flip the "
"image along the first axis before running Crop",
img->name, p->cp.hdu);
/* 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->pixscale)
{
for(i=0;ipixscale[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);
}
else
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 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;isized[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'.
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 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 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", 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],
img->corners[6], img->corners[7]);
exit(0);
*/
}
/*******************************************************************/
/************ Check if WCS is in image **************/
/*******************************************************************/
/* Set the four sides around the point of interest in RA and Dec.
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
wcsmode_crop_corners(struct onecropparams *crp)
{
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)
{
/* Find their minimum and maximum values. */
for(i=0;invertices;++i)
{
if(p->wpolygon[i*2]>maxra) maxra=p->wpolygon[i*2];
if(p->wpolygon[i*2]wpolygon[i*2];
if(p->wpolygon[i*2+1]>maxdec) maxdec=p->wpolygon[i*2+1];
if(p->wpolygon[i*2+1]wpolygon[i*2+1];
}
/* Set the corners: */
crp->corners[0] = maxra; crp->corners[1] = mindec; /* Bottom Left */
crp->corners[2] = minra; crp->corners[3] = mindec; /* Bottom Right */
crp->corners[4] = maxra; crp->corners[5] = maxdec; /* Top Left */
crp->corners[6] = minra; crp->corners[7] = maxdec; /* Top Right */
}
else
{
/* 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];
/* Calculate the declination in radians for easy readability. */
dr=d*M_PI/180;
/* 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 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 above `point_in_dataset'. */
if(crp->corners[1]*(crp->corners[1]+crp->sized[1]) < 0 )
{
/* 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 above `point_in_dataset'. */
crp->equatorcorr[1]=crp->sized[0]*cos(crp->corners[1]*M_PI/180);
}
/* Just to check:
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);
*/
}
/* We have the polygon vertices in WCS coordinates and need to change them
to one input image's pixel coordinates. */
void
fillcrpipolygon(struct onecropparams *crp)
{
struct cropparams *p=crp->p;
gal_data_t *tmp, *coords=NULL;
size_t i, d, ndim=p->imgs->ndim;
/* Allocate the necessary arrays for each column. */
for(d=0;dnvertices,
NULL, 0, -1, NULL, NULL, NULL);
/* Fill in the world coordinate columns. */
for(i=0;invertices;++i)
{
d=0;
for(tmp=coords;tmp!=NULL;tmp=tmp->next)
((double *)(tmp->array))[i] = p->wpolygon[ i * ndim + d++ ];
}
/* Convert them to image coordinates. */
gal_wcs_world_to_img(coords, p->imgs[crp->in_ind].wcs, 1);
/* Allocate the image polygon array, and put the image polygon vertice
values into it. */
crp->ipolygon=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim*p->nvertices, 0,
__func__, "crp->ipolygon");
for(i=0;invertices;++i)
{
d=0;
for(tmp=coords;tmp!=NULL;tmp=tmp->next)
crp->ipolygon[ i * ndim + d++ ] = ((double *)(tmp->array))[i];
}
/* Clean up. */
gal_list_data_free(coords);
}
/* BASICS:
=======
An image is a rectangle, but the sky is on a globe. When the images
are aligned to the celstial coordinates, (as we have required here
in wcscheckprepare) the first FITS axis shows change in RA, while
the second axis shows change in Dec. The declination always changes
along a great circle, so there is no problem. But unless
declination is constrained to zero, RA changes on small circles.
See the rectangle below, assume it is an image. To check if a given
point is within the same declination as this rectangle is very
simple, since d3==d4 and d1==d2. Note that r1>r2 and r3>r4 (because
right ascension increases to the east).
(r3,d3) ------------------ (r4,d4)
| |
| |
| |
| |
(r1,d1) ------------------ (r2,d2)
But for RA, the same number of pixels on each declination,
corresponds to different ranges in right ascention. As the
declination gets higher in the northern hemisphere (where the
declination rises towards the top of the image) r1-r2 becomes
smaller than r3-r4. So, in terms of coverage in RA and Dec, this
box should rather be shown like this trapezoid (exaggerated):
--------------------
| |
| | (Northern hemisphere)
| |
| |
--------------
On the southern hemisphere it should be shown like this:
----------------
| |
| | (Southern hemisphere)
| |
| |
----------------------
The functional form of the change is the inverse cosine, so:
(r3-r4)=(r1-r2)/cos(d3-d1) (North)
(r1-r2)=(r3-r4)/cos(d1-d3) (South)
QUESTION:
========
Is a given point at the RA and Dec of (rp,dp) inside our
rectangular image?
IMAGE IS FULLY WITHIN ONE HEMISPHERE
------------------------------------
Our reference point for the image is the first pixel in the image,
which by convention is the (r1,d1) point in the rectangle above. We
also have the angular size of the rectangular image as 'sr', 'sd'
(for "size in RA" and "size in dec"). NOTE: This has to be We also
assume r1+sr and d1+sd are the distances to the last pixels in our
rectangular image.
As explained above, to check the declination range, everything is
very easy:
dp>=d1 && dp<=d1+sd
For RA, things become a little more complicated (recall that
r1>r3). `n` is defined as half of the extra space between the top
and bottom lines of the two trapezoids.
(North) n=0.5*sr*(1/cos(dp-d1)-1) ==> rp<=r1+n && rp>=r1-sr-n
(South) n=0.5*sr*(1-cos(dp-d1)) ==> rp<=r1-n && rp>=r1-sr+n
IMAGE CROSSES THE EQUATOR
-------------------------
When d1*(d1+sd)<0, the image crosses the equator (d1 is negative
and d1+sd is positive). In this case, we define `re` and `sre` as
an equivalent of r1 and sr but on the equator:
re=r1-0.5*sr*(1-cos(d1)) && sre=sr*cos(d1)
then for all the points with negative declination we use the
(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.
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;
/* 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! */
n=0.5f*s[0]*( 1 - cos((p[1]-i[1])*M_PI/180) );
if(p[0]<=i[0]-n && p[0]>=i[0]-s[0]+n)
return 1;
}
else /* Point is in the northern hemisphere. */
{
if( i[1] * (s[1]+i[1]) > 0 ) /* Image does not cross equator. */
{
n=0.5f*s[0]*( 1/cos((p[1]-i[1])*M_PI/180) - 1);
if(p[0]<=i[0]+n && p[0]>=i[0]-s[0]-n)
return 1;
}
else /* Image crosses the equator.*/
{
n=0.5f*c[1]*( 1/cos((p[1] )*M_PI/180) - 1);
if(p[0]<=c[0]+n && p[0]>=c[0]-c[1]-n)
return 1;
}
}
}
return 0;
}
/* Is there an overlap between this crop box and the survey image?
This function will return 0 if there isn't and 1 if there is.
We don't want to necessarily assume that the crop box is smaller
than the survey images. If we made that assumption, we only had to
check if the corners of the crop are in the image. When we allow
the input survey images to be smaller than the crop box, it might
happen that none of the corners of the crop are in the image but
there is an overlap (the survey image is completely within the crop
box). So we have to check both. */
int
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 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
{
/* 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(dsized;
i=crp->corners;
c=crp->equatorcorr;
fd=(d=p->imgs[crp->in_ind].corners) + 8;
do
{
if( point_in_dataset(d, i, s, c, ndim) ) return 1;
d+=ndim;
}
while(d