[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master c3bd7956: Library (polygon.h): new function to
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master c3bd7956: Library (polygon.h): new function to find polygon area on sky |
|
Date: |
Sat, 2 Sep 2023 19:20:28 -0400 (EDT) |
branch: master
commit c3bd7956f1f458bed3300c2384535b428a0fcf91
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (polygon.h): new function to find polygon area on sky
Until now, the only available function to find the area of a pixel was the
'gal_polygon_area' function, which assumed a flat coordinate
system. However, this function was used in the 'gal_warp_pixelarea'
function (that calculates the pixel are on the spherical sky
coordinates!). As a result, the reported areas were wrong (would change
with declination).
With this commit, that function has been renamed to 'gal_polygon_area_flat'
to remind the callers that it is for a flat coordinate system. For
spherical coordinate system, a new 'gal_polygon_area_sky' has been
added.
Also, since the area of the pixels are independent of the position on the
sky (CRVALi in the FITS keywords), the 'gal_warp_pixelarea' now does the
WCS conversion with a temporary CRVALi of (180,0). To do this, it was
necessary to define a new function to change CRVAL (because WCSLIB keeps
and uses low-level internal parmaters that will not change when we simply
change 'wcs->crval').
This fixes bug #64615.
---
NEWS | 8 ++++
bin/fits/meta.c | 2 +-
bin/warp/warp.c | 7 +--
doc/gnuastro.texi | 31 +++++++++++++-
lib/gnuastro/polygon.h | 5 ++-
lib/gnuastro/wcs.h | 3 ++
lib/polygon.c | 75 +++++++++++++++++++++++++-------
lib/warp.c | 114 ++++++++++++++++++++++++++++++++++---------------
lib/wcs.c | 59 ++++++++++++++++++++++---
9 files changed, 241 insertions(+), 63 deletions(-)
diff --git a/NEWS b/NEWS
index 9c95b923..43c09116 100644
--- a/NEWS
+++ b/NEWS
@@ -81,6 +81,7 @@ See the end of the file for license conditions.
- gal_blank_not_minmax_coords: returns the minimum/maximum coordinates of
non-blank regions of input dataset.
- gal_blank_trim: trim all outer blank regions from the input dataset.
+ - gal_polygon_area_sky: area of polygon defined on celestial coordinates.
- gal_pool_min: min-pooling function, see 'pool-min' above.
- gal_pool_max: max-pooling function, see 'pool-min' above.
- gal_pool_sum: sum-pooling function, see 'pool-min' above.
@@ -157,6 +158,12 @@ See the end of the file for license conditions.
would cause problem in the final stack and were not removed until now.
Implemented by Sepideh Eskandarlou
+ Library:
+ - gal_polygon_area_flat: new name for the old 'gal_polygon_area'
+ function. This was necessary because this function assumes flat
+ coordinates and can lead to wrong results when used with a polygon
+ defined by RA,Dec for example.
+
** Bugs fixed
bug #64138: Arithmetic's mknoise-poisson only using first pixel value.
@@ -195,6 +202,7 @@ See the end of the file for license conditions.
doesn't exist. Reported by Sepideh Eskandarlou.
bug #64610: --copykeys of Fits program finds wrong key when the first
characters of the keyword names match.
+ bug #64615: Fits program's '--pixelareaonwcs' only accurate on equator.
diff --git a/bin/fits/meta.c b/bin/fits/meta.c
index 64fc526e..d1848657 100644
--- a/bin/fits/meta.c
+++ b/bin/fits/meta.c
@@ -148,7 +148,7 @@ meta_pixelareaonwcs(struct fitsparams *p)
/* Call for an empty wcsalign structure. */
gal_warp_wcsalign_t wa=gal_warp_wcsalign_template();
- /* Initialize the warping variables based on commandline arguments. */
+ /* Initialize the warping variables based on command-line arguments. */
meta_initialize(p, &wa);
/* Execute the warping and fill the data-structure with results. */
diff --git a/bin/warp/warp.c b/bin/warp/warp.c
index d8cad9a5..4d025e5e 100644
--- a/bin/warp/warp.c
+++ b/bin/warp/warp.c
@@ -148,7 +148,8 @@ warp_onthread_linear(void *inparam)
ocrn[j*2], ocrn[j*2+1], icrn_base[j*2],
icrn_base[j*2+1]);
printf("------- Ordered -------\n");
- for(j=0;j<4;++j) printf("(%.3f, %.3f)\n", icrn[j*2], icrn[j*2+1]);
+ for(j=0;j<4;++j)
+ printf("(%.3f, %.3f)\n", icrn[j*2], icrn[j*2+1]);
printf("------- Start and ending pixels -------\n");
printf("X: %ld -- %ld\n", xstart, xend);
printf("Y: %ld -- %ld\n", ystart, yend);
@@ -177,7 +178,7 @@ warp_onthread_linear(void *inparam)
/* Find the overlapping (clipped) polygon: */
gal_polygon_clip(icrn, 4, pcrn, 4, ccrn, &numcrn);
- area=gal_polygon_area(ccrn, numcrn);
+ area=gal_polygon_area_flat(ccrn, numcrn);
/* Add the fractional value of this pixel. If this
output pixel covers a NaN pixel in the input grid,
@@ -367,7 +368,7 @@ warp_linear_init(struct warpparams *p)
forarea[2*i]=icrn[2*p->ordinds[i]];
forarea[2*i+1]=icrn[2*p->ordinds[i]+1];
}
- p->opixarea=gal_polygon_area(forarea, 4);
+ p->opixarea=gal_polygon_area_flat(forarea, 4);
/* Find which index after transformation will have the minimum and
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index e4038118..6e44a27c 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -38181,6 +38181,12 @@ If you just want to write the WCS into an empty file,
you can use @code{gal_wcs_
Return a fully allocated (independent) copy of @code{wcs}.
@end deftypefun
+@deftypefun {struct wcsprm *} gal_wcs_copy_new_crval (struct wcsprm
@code{*wcs}, double @code{*crval})
+Return a fully allocated (independent) copy of @code{wcs} with a new set of
@code{CRVAL} values.
+WCSLIB keeps a lot of extra information within @code{wcsprm} and for
optimizations, those extra information are used in its calculations.
+Therefore, if you want to change parameters like the reference point's sky
coordinate values (@code{CRVAL}), simply changing the values in
@code{wcs->crval[0]} or @code{wcs->crval[1]} will not affect WCSLIB's
calculations; you need to call this function.
+@end deftypefun
+
@deftypefun void gal_wcs_remove_dimension (struct wcsprm @code{*wcs}, size_t
@code{fitsdim})
Remove the given FITS dimension from the given @code{wcs} structure.
@end deftypefun
@@ -39465,9 +39471,30 @@ Returns @code{1} if the polygon is convex with
vertices defined by @code{v} and
Note that the vertices of the polygon should be sorted in an anti-clockwise
manner.
@end deftypefun
-@deftypefun double gal_polygon_area (double @code{*v}, size_t @code{n})
-Find the area of a polygon with vertices defined in @code{v}.
+@deftypefun double gal_polygon_area_flat (double @code{*v}, size_t @code{n})
+Find the area of a polygon with vertices defined in @code{v} on a euclidian
(flat) coordinate system.
+@code{v} points to an array of doubles which keep the positions of the
vertices such that @code{v[0]} and @code{v[1]} are the positions of the first
vertice to be considered.
+@end deftypefun
+
+@deftypefun double gal_polygon_area_sky (double @code{*v}, size_t @code{n})
+Find the area of a polygon with vertices defined in @code{v} on a celestial
coordinate system.
+This is a coordinate system where the first coordinate goes from 0 to 360
(increasing to the right), while the second coordinate ranges from -90 to +90
(on the poles).
@code{v} points to an array of doubles which keep the positions of the
vertices such that @code{v[0]} and @code{v[1]} are the positions of the first
vertice to be considered.
+
+This function uses an approximation to account for the curvature of the sky
and the different nature of spherical coordinates with respect to the flat
coordinate system.
+@url{https://savannah.gnu.org/bugs/index.php?64617, Bug 64617} has been
defined in Gnuastro to address this problem.
+Please check that bug in case it has been fixed.
+Until this bug is fixed, here are some tips:
+@itemize
+@item
+Subtract the RA and Dec of all the vertice coordinates from a constant so the
center of the polygon falls on (RA, Dec) of (180,0).
+The sphere has a similar nature everywhere on it, so shifting the polygon
vertices will not change its area; this also removes issues with the RA=0 or
RA=360 coordinate and decrease issues caused by RA depending on declination.
+@item
+These approximations should not cause any statistically significant error on
normal (less than a few degrees) scales.
+But it won't hurt to do a small sanity check for your particular usage
scenario.
+@item
+Any help (even in the mathematics of the problem; not necessary programming)
would be appreciated (we didn't have time to derive the necessary equations),
so if you have some background in this and can prepare the mathematical
description of the problem, please get in touch.
+@end itemize
@end deftypefun
@deftypefun int gal_polygon_is_inside (double @code{*v}, double @code{*p},
size_t @code{n})
diff --git a/lib/gnuastro/polygon.h b/lib/gnuastro/polygon.h
index 2dfef7d9..9f891d03 100644
--- a/lib/gnuastro/polygon.h
+++ b/lib/gnuastro/polygon.h
@@ -65,7 +65,10 @@ int
gal_polygon_is_convex(double *v, size_t n);
double
-gal_polygon_area(double *v, size_t n);
+gal_polygon_area_flat(double *v, size_t n);
+
+double
+gal_polygon_area_sky(double *v, size_t n);
int
gal_polygon_is_inside(double *v, double *p, size_t n);
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 5e46b5a6..17a18d66 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -177,6 +177,9 @@ gal_wcs_distortion_convert(struct wcsprm *inwcs,
struct wcsprm *
gal_wcs_copy(struct wcsprm *wcs);
+struct wcsprm *
+gal_wcs_copy_new_crval(struct wcsprm *in, double *crval);
+
void
gal_wcs_remove_dimension(struct wcsprm *wcs, size_t fitsdim);
diff --git a/lib/polygon.c b/lib/polygon.c
index 91540cb1..479d9b41 100644
--- a/lib/polygon.c
+++ b/lib/polygon.c
@@ -47,14 +47,18 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* The cross product of two points from the center. */
#define GAL_POLYGON_CROSS_PRODUCT(A, B) ( (A)[0]*(B)[1] - (B)[0]*(A)[1] )
-
-
+/* For spherical coordinates: the RA is only in units of degrees on the
+ equator: as we got to higher or lower declinations, any difference in RA
+ should be multiplied by the cos(DEC).*/
+#define GAL_POLYGON_CROSS_PRODUCT_SKY(A, B) \
+ ( (A)[0]*(B)[1]*cos((A)[1]*M_PI/180) \
+ - (B)[0]*(A)[1]*cos((B)[1]*M_PI/180) )
/* Find the cross product (2*area) between three points. Each point is
assumed to be a pointer that has atleast two values within it. */
-#define GAL_POLYGON_TRI_CROSS_PRODUCT(A, B, C) \
- ( ( (B)[0]-(A)[0] ) * ( (C)[1]-(A)[1] ) - \
- ( (C)[0]-(A)[0] ) * ( (B)[1]-(A)[1] ) ) \
+#define GAL_POLYGON_TRI_CROSS_PRODUCT(A, B, C) \
+ ( ( (B)[0]-(A)[0] ) * ( (C)[1]-(A)[1] ) \
+ - ( (C)[0]-(A)[0] ) * ( (B)[1]-(A)[1] ) )
@@ -130,11 +134,12 @@ gal_polygon_vertices_sort_convex(double *in, size_t n,
size_t *ordinds)
tindexs[GAL_POLYGON_MAX_CORNERS];
if(n>GAL_POLYGON_MAX_CORNERS)
- error(EXIT_FAILURE, 0, "%s: most probably a bug! The number of corners "
- "is more than %d. This is an internal value and cannot be set from "
- "the outside. Most probably some bug has caused this un-normal "
- "value. Please contact us at %s so we can solve this problem",
- __func__, GAL_POLYGON_MAX_CORNERS, PACKAGE_BUGREPORT);
+ error(EXIT_FAILURE, 0, "%s: most probably a bug! The number of "
+ "corners is more than %d. This is an internal value and "
+ "cannot be set from the outside. Most probably some bug "
+ "has caused this un-normal value. Please contact us at %s "
+ "so we can solve this problem", __func__,
+ GAL_POLYGON_MAX_CORNERS, PACKAGE_BUGREPORT);
/* For a check:
printf("\n\nBefore sorting:\n");
@@ -234,12 +239,11 @@ gal_polygon_is_convex(double *v, size_t n)
vertices such that v[0] and v[1] are the positions of the first
corner to be considered.
- We will start from the edge connecting the last pixel to the first
- pixel for the first step of the loop, for the rest, j is always
- going to be one less than i.
- */
+ To optimize the process (avoid repetition): We start from the edge
+ connecting the last pixel to the first pixel for the first step of the
+ loop, for the rest, j is always going to be one less than i. */
double
-gal_polygon_area(double *v, size_t n)
+gal_polygon_area_flat(double *v, size_t n)
{
double sum=0.0f;
size_t i=0, j=n-1;
@@ -256,6 +260,47 @@ gal_polygon_area(double *v, size_t n)
+/* Spherical coordinates need special attention.*/
+double
+gal_polygon_area_sky(double *v, size_t n)
+{
+ int a1b0;
+ size_t i=0, j=n-1;
+ double *a, *b, sum=0.0f;
+
+ /* Go over all the vertices of the polygon. */
+ while(i<n)
+ {
+ /* For easy reading. */
+ a=v+j*2;
+ b=v+i*2;
+
+ /* We need to account for passing the two points pass over the line
+ of RA=0 (assuming that the vertices are not larger than 180
+ degrees apart!) */
+ if( fabs(a[0]-b[0]) < 180.0f)
+ sum+=GAL_POLYGON_CROSS_PRODUCT_SKY(a, b);
+ else
+ {
+ /* Add the smaller RA by 360. */
+ if(a[0]<b[0]) {a1b0=1; a[0]+=360.0f;}
+ else {a1b0=0; b[0]+=360.0f;}
+
+ /* Add the cross product to the sum. */
+ sum+=GAL_POLYGON_CROSS_PRODUCT_SKY(a, b);
+
+ /* Return the changed value. */
+ if(a1b0) a[0]-=360.0f; else b[0]-=360.0f;
+ }
+ j=i++;
+ }
+ return fabs(sum)/2.0f;
+}
+
+
+
+
+
/* This fuction test if a point is inside the polygon using winding number
algorithm. See its wiki here:
https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm
diff --git a/lib/warp.c b/lib/warp.c
index 55e8e2b7..621d2b2a 100644
--- a/lib/warp.c
+++ b/lib/warp.c
@@ -77,7 +77,6 @@ warp_alloc_perimeter(gal_data_t *input)
/* High level variables */
size_t npcrn=2*(is0+is1);
-
pcrn=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &npcrn, NULL, 0,
minmapsize, quietmmap, NULL, NULL, NULL);
pcrn->next=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &npcrn, NULL, 0,
@@ -117,7 +116,7 @@ warp_alloc_perimeter(gal_data_t *input)
/* Create a base image with WCS consisting of the basic geometry
- keywords */
+ keywords. */
static void
warp_wcsalign_init_output_from_params(gal_warp_wcsalign_t *wa)
{
@@ -257,13 +256,14 @@ warp_wcsalign_init_output_from_params(gal_warp_wcsalign_t
*wa)
+/* Initialize the vertices of each output pixel (accounting for any
+ edge-sampling). */
static void
warp_wcsalign_init_vertices(gal_warp_wcsalign_t *wa)
{
size_t es=wa->edgesampling;
size_t ind, i, j, ix, iy;
- gal_data_t *vertices=NULL;
double gap=1.0f/(es+1.0f);
size_t os0=wa->output->dsize[0];
size_t os1=wa->output->dsize[1];
@@ -276,22 +276,26 @@ warp_wcsalign_init_vertices(gal_warp_wcsalign_t *wa)
size_t v0=WARP_WCSALIGN_V0(es,os0,os1);
size_t nvertices=nvcrn*(os1+1)+nhcrn*(os0+1);
- /* Now create all subpixels based on the edgesampling '-E' option */
- vertices=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nvertices, NULL, 0,
- minmapsize, quietmmap, "OutputRA", NULL, NULL);
- vertices->next=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nvertices,
- NULL, 0, minmapsize, quietmmap,
- "OutputDec", NULL, NULL);
- x=vertices->array; y=vertices->next->array;
-
+ /* Allocate the space for all the vertice coordinates. */
+ wa->vertices=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nvertices,
+ NULL, 0, minmapsize, quietmmap, NULL,
+ NULL, NULL);
+ wa->vertices->next=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nvertices,
+ NULL, 0, minmapsize, quietmmap,
+ NULL, NULL, NULL);
+
+ /* Parse all pixels. */
+ x=wa->vertices->array;
+ y=wa->vertices->next->array;
for(ind=os0*os1; ind--;)
{
+ /* For easy reading. */
row=ind%os1;
col=floor(ind/os1);
ix=WARP_WCSALIGN_H(ind,es,os1);
iy=WARP_WCSALIGN_V(ind,es,v0,os1);
- /* Bottom left */
+ /* Bottom left edge of pixel. */
x[ix]= 0.5f+row;
y[ix]= 0.5f+col;
@@ -320,10 +324,11 @@ warp_wcsalign_init_vertices(gal_warp_wcsalign_t *wa)
/* Right */
for(ind=os1-1; ind<os1*os0; ind+=os1)
{
- row=(size_t)(ind%os1); col=(size_t)(ind/os1);
-
- iy=WARP_WCSALIGN_V(ind,es,v0,os1);
+ /* For easy reading. */
+ row=(size_t)(ind%os1);
+ col=(size_t)(ind/os1);
ix=WARP_WCSALIGN_H(ind,es,os1);
+ iy=WARP_WCSALIGN_V(ind,es,v0,os1);
/* Bottom right */
j=ix+es+1;
@@ -339,22 +344,28 @@ warp_wcsalign_init_vertices(gal_warp_wcsalign_t *wa)
}
}
- wa->vertices=vertices;
+ /* For a check:
+ for(i=0;i<nvertices;++i)
+ printf("%zu: %f, %f\n", i, x[i], y[i]);
+ printf("%s: GOOD\n", __func__); exit(0);
+ */
}
+/* Check if the vertice orientation is clock-wise or counter-clock-wise. */
static void
-warp_check_output_orientation(gal_warp_wcsalign_t *wa)
+warp_check_output_clockwise(gal_warp_wcsalign_t *wa)
{
size_t gcrn=wa->gcrn;
size_t es=wa->edgesampling;
double *vx=wa->vertices->array;
double *vy=wa->vertices->next->array;
size_t indices[4]={ 0, es+1, gcrn+es+1, gcrn };
- double *temp=gal_pointer_allocate(GAL_TYPE_FLOAT64, 8, 0, __func__, NULL);
+ double *temp=gal_pointer_allocate(GAL_TYPE_FLOAT64, 8, 0, __func__,
+ NULL);
temp[ 0 ]=vx[ indices[0] ]; temp[ 1 ]=vy[ indices[0] ];
temp[ 2 ]=vx[ indices[1] ]; temp[ 3 ]=vy[ indices[1] ];
@@ -585,8 +596,7 @@ warp_wcsalign_init_output_from_wcs(gal_warp_wcsalign_t *wa,
-/* Check parameters shared between the wcsalign procedure and the other
- library functions such as pixelarea_onwcs. */
+/* Check parameters for aligning and pixelarea. */
static void
warp_check_basic_params(gal_warp_wcsalign_t *wa, const char *func)
{
@@ -633,13 +643,11 @@ warp_check_basic_params(gal_warp_wcsalign_t *wa, const
char *func)
static void
warp_wcsalign_init_internals(gal_warp_wcsalign_t *wa)
{
- size_t es, os0, os1;
-
- es=wa->edgesampling;
- os0=wa->output->dsize[0];
- os1=wa->output->dsize[1];
+ size_t es=wa->edgesampling;
+ size_t os0=wa->output->dsize[0];
+ size_t os1=wa->output->dsize[1];
- /* Pickle variables so other functions can access them. */
+ /* Common variables to simplify next functions. */
wa->ncrn=4*es+4;
wa->gcrn=1+os1*(es+1);
wa->v0=WARP_WCSALIGN_V0(es, os0, os1);
@@ -647,7 +655,7 @@ warp_wcsalign_init_internals(gal_warp_wcsalign_t *wa)
/* Determine the output image rotation direction so we can sort the
indices in counter clockwise order. This is necessary for the
'gal_polygon_clip' function to work. */
- warp_check_output_orientation(wa);
+ warp_check_output_clockwise(wa);
}
@@ -967,10 +975,18 @@ gal_warp_wcsalign_onpix(gal_warp_wcsalign_t *wa, size_t
ind)
/* Read the value of the input pixel. */
v=inputarr[(y-1)*is1+x-1];
- /* Find the overlapping (clipped) polygon: */
+ /* Find the overlapping (clipped) polygon and its area.
+
+ In theory, instead of 'gal_polygon_area_flat', we should be
+ using the spherical polygon area ('gal_polygon_area_flat').
+ But the area that is calculate here is not used in an absolute
+ sense: it is only relative for comparison with other input
+ pixels that overlap with this output pixel. Therefore, because
+ the flat area calculation is faster, we'll suffice to that
+ unless we discover there is any problem with it. */
numcrn=0; /* initialize it. */
gal_polygon_clip(ocrn, ncrn, pcrn, 4, ccrn, &numcrn);
- area=gal_polygon_area(ccrn, numcrn);
+ area=gal_polygon_area_flat(ccrn, numcrn);
/* Write each pixel's maximum coverage fraction if asked. */
if( maxfrac ) maxfrac[ind] = fmax(area, maxfrac[ind]);
@@ -998,9 +1014,11 @@ gal_warp_wcsalign_onpix(gal_warp_wcsalign_t *wa, size_t
ind)
if( maxfrac && maxfrac[ind]==-DBL_MAX ) maxfrac[ind]=NAN;
/* See if the pixel value should be set to NaN or not (because of not
- enough coverage). Note that 'ocrn' is sorted in anti-clockwise
- order already. */
- opixarea=gal_polygon_area(ocrn, ncrn);
+ enough coverage). Note that 'ocrn' is sorted in anti-clockwise order
+ already. For a description of why we are not using
+ 'gal_polygon_area_sky', see the comment above the previous call to
+ 'gal_polygon_area_flat' above. */
+ opixarea=gal_polygon_area_flat(ocrn, ncrn);
if( numinput && filledarea/opixarea < wa->coveredfrac-1e-5)
numinput=0;
@@ -1134,13 +1152,14 @@ warp_pixelarea_onthread(void *inparam)
/* Loop over pixels given from the 'warp' function */
for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
{
+ /* Pixel to use. */
ind=tprm->indexs[i];
/* Fix the vertice ordering, crucial for calculating the area. */
ocrn=warp_pixel_perimeter(wa, ind);
/* Now that the vertices are in CCW order, calculate the area. */
- area=gal_polygon_area(ocrn, wa->ncrn);
+ area=gal_polygon_area_sky(ocrn, wa->ncrn);
outputarr[ind]=area;
/* Clean up the allocated array. */
@@ -1162,8 +1181,11 @@ warp_pixelarea_onthread(void *inparam)
void
gal_warp_pixelarea(gal_warp_wcsalign_t *wa)
{
+ struct wcsprm *wcs;
gal_data_t *input=wa->input;
+ double crval[2]={180.0f,0.0f};
+ /* Basic sanity check. */
warp_check_basic_params(wa, __func__);
/* Create the output dataset. */
@@ -1175,7 +1197,31 @@ gal_warp_pixelarea(gal_warp_wcsalign_t *wa)
/* Create the vertices based on the edgesampling value. */
warp_wcsalign_init_vertices(wa);
warp_wcsalign_init_internals(wa);
- gal_wcs_img_to_world(wa->vertices, input->wcs, 1);
+
+ /* For a check (before conversion).
+ size_t i;
+ double *x=wa->vertices->array;
+ double *y=wa->vertices->next->array;
+ for(i=0;i<wa->vertices->size;++i)
+ printf("%-20.15f %-20.15f\n", x[i], y[i]);
+ */
+
+ /* Change the CRVALs to be on 0 and 180 (this helps to avoid issues with
+ an image passing the RA=0.0, or high declination problems. */
+ wcs=gal_wcs_copy_new_crval(input->wcs, crval);
+
+ /* Convert the output pixel-coordinate vertices to WCS and free the
+ temporary WCS. */
+ gal_wcs_img_to_world(wa->vertices, wcs, 1);
+ gal_wcs_free(wcs);
+
+ /* For a check (after conversion).
+ size_t j;
+ double *ra=wa->vertices->array;
+ double *dec=wa->vertices->next->array;
+ for(j=0;j<wa->vertices->size;++j)
+ printf("%-20.15f %-20.15f\n", ra[j], dec[j]);
+ */
/* Calculate pixel area on WCS and write to output. */
gal_threads_spin_off(warp_pixelarea_onthread, wa, wa->output->size,
diff --git a/lib/wcs.c b/lib/wcs.c
index 1c8fc697..86c52ecf 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -1255,8 +1255,8 @@ gal_wcs_distortion_convert(struct wcsprm *inwcs, int
outdisptype,
outwcs=gal_wcsdistortion_sip_to_tpv(inwcs);
break;
default:
- error(EXIT_FAILURE, 0, "%s: conversion from %s to %s is not yet "
- "supported. Please contact us at '%s'", __func__,
+ error(EXIT_FAILURE, 0, "%s: conversion from %s to %s is not "
+ "yet supported. Please contact us at '%s'", __func__,
gal_wcs_distortion_to_string(indisptype),
gal_wcs_distortion_to_string(outdisptype),
PACKAGE_BUGREPORT);
@@ -1274,8 +1274,8 @@ gal_wcs_distortion_convert(struct wcsprm *inwcs, int
outdisptype,
outwcs=gal_wcsdistortion_tpv_to_sip(inwcs, fitsize);
break;
default:
- error(EXIT_FAILURE, 0, "%s: conversion from %s to %s is not yet "
- "supported. Please contact us at '%s'", __func__,
+ error(EXIT_FAILURE, 0, "%s: conversion from %s to %s is not "
+ "yet supported. Please contact us at '%s'", __func__,
gal_wcs_distortion_to_string(indisptype),
gal_wcs_distortion_to_string(outdisptype),
PACKAGE_BUGREPORT);
@@ -1293,9 +1293,9 @@ gal_wcs_distortion_convert(struct wcsprm *inwcs, int
outdisptype,
/* A bug! This distortion is not yet recognized. */
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
- "the problem. The identifier '%d' is not recognized as a "
- "distortion", __func__, PACKAGE_BUGREPORT, indisptype);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The identifier '%d' is not recognized "
+ "as a distortion", __func__, PACKAGE_BUGREPORT, indisptype);
}
/* Return the converted WCS. */
@@ -1370,6 +1370,51 @@ gal_wcs_copy(struct wcsprm *wcs)
+/* Copy the given WCS into a new one with differnet CRVALs. Because WCSLIB
+ keeps internal constructs to speed up operations, we cannot simply
+ change these values, we need to write the WCS into a set of keywords,
+ then read them as a new one. */
+struct wcsprm *
+gal_wcs_copy_new_crval(struct wcsprm *in, double *crval)
+{
+ char *wcsstr;
+ double *ocrval;
+ int status=0, relax=WCSHDR_all;
+ struct wcsprm *incpy, *out=NULL;
+ int nkeys, nwcs, ctrl=0, nreject=0;
+
+ /* Copy the input WCS structure into a new one (in case the caller is
+ using it in another function). */
+ incpy=gal_wcs_copy(in);
+
+ /* Keep the original pointer of CRVAL: for some reason 'wcsprm' cannot
+ deal with a 'NULL' CRVAL, so we need to return the original CRVAL to
+ the 'wcsprm' before freeing it (we do not want to free the caller's
+ CRVAL). */
+ ocrval=incpy->crval;
+ incpy->crval=crval;
+
+ /* Convert the WCS into a string and read it into a new WCS structure. We
+ do not need all the checks in 'gal_wcs_read_fitsptr' because this
+ string was written just now. */
+ wcsstr=gal_wcs_write_wcsstr(incpy, &nkeys);
+ status=wcspih(wcsstr, nkeys, relax, ctrl, &nreject, &nwcs, &out);
+ if(status)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to fix the "
+ "problem. The internally created WCS string could not be parsed "
+ "as a new WCS structure", __func__, PACKAGE_BUGREPORT);
+
+ /* Clean up and return. */
+ incpy->crval=ocrval;
+ gal_wcs_free(incpy);
+ free(wcsstr);
+ return out;
+}
+
+
+
+
+
/* Remove the algorithm part of CTYPE (anything after, and including, a
'-') if necessary. */
static void
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master c3bd7956: Library (polygon.h): new function to find polygon area on sky,
Mohammad Akhlaghi <=