[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 9321009: Decomposing PCi_j and CDELTi matices
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 9321009: Decomposing PCi_j and CDELTi matices in output WCS |
Date: |
Wed, 18 Jan 2017 01:55:44 +0000 (UTC) |
branch: master
commit 93210093628f65a271190b46d5d72a6d078a5202
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Decomposing PCi_j and CDELTi matices in output WCS
When reading images that have their WCS format in the `CDi_j' formalism,
WCSLIB will simply put the `CDi_j' array into the `PCi_j' array and set all
the `CDELTi' elements to `1.0'.
According to the FITS standard, in the `PCi_j' formalism, the scale factors
are encoded in the `CDELTi' keywords and the `PCi_j' matrix should contain
all the other transformations (recall that scaling is just one of all the
possible linear transformations).
Therefore, as noted by Lee Kelvin in bug #50073, this behavior of WCSLIB
can cause problems for programs that look into the value of `CDELTi' array
independent of the `PCi_j' matrix (like ds9 when viewing two files).
Fortunately, as part of the solution to bug #50072 (which was also reported
by Lee), a robust way to find the pixel scale was found using `Singular
Value Decomposition'. So when the `wcsprm' structure has `CDELTi' values of
`1.0', the new `gal_wcs_decompose_pc_cdelt' function is called to decompose
the two arrays so WCSLIB will print realistic `CDELTi' keywords.
These corrections also allowed for more clear pixel resolution checks in
ImageCrop's WCS mode crops.
This fixes bug #50073.
---
bin/imgcrop/wcsmode.c | 22 +++++++++++++-------
lib/fits.c | 4 ++++
lib/gnuastro/wcs.h | 3 +++
lib/wcs.c | 55 +++++++++++++++++++++++++++++++++++++++++++++++++
4 files changed, 76 insertions(+), 8 deletions(-)
diff --git a/bin/imgcrop/wcsmode.c b/bin/imgcrop/wcsmode.c
index ae894d8..0072d9a 100644
--- a/bin/imgcrop/wcsmode.c
+++ b/bin/imgcrop/wcsmode.c
@@ -50,7 +50,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
void
wcscheckprepare(struct imgcropparams *p, struct inputimgs *img)
{
- double twidth;
+ double twidth, *pixscale;
struct wcsprm *wcs=img->wcs;
int i, status[4]={0,0,0,0}, ncoord=4, nelem=2;
double imgcrd[8], phi[4], theta[4], pixcrd[8];
@@ -88,7 +88,9 @@ wcscheckprepare(struct imgcropparams *p, struct inputimgs
*img)
if(p->res==0.0f)
{
/* Set the resolution of the image. */
- p->res=wcs->pc[3];
+ pixscale=gal_wcs_pixel_scale_deg(wcs);
+ p->res=pixscale[0];
+ free(pixscale);
/* Set the widths such that iwidth and wwidth are exactly the same
(within their different units ofcourse). Also make sure that the
@@ -109,12 +111,16 @@ wcscheckprepare(struct imgcropparams *p, struct inputimgs
*img)
p->iwidth[1]=p->iwidth[0];
}
else
- if(p->res!=wcs->pc[3])
- 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);
+ {
+ pixscale=gal_wcs_pixel_scale_deg(wcs);
+ 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);
+ }
/* Get the coordinates of the first pixel in the image. */
diff --git a/lib/fits.c b/lib/fits.c
index be47af8..429f026 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <unistd.h>
#include <gnuastro/git.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
#include "checkset.h"
@@ -1848,6 +1849,9 @@ gal_fits_array_to_file(char *filename, char *extname, int
bitpix, void *array,
if(wcs)
{
+ /* Decompose the `PCi_j' matrix and `CDELTi' vector. */
+ gal_wcs_decompose_pc_cdelt(wcs);
+
/* Convert the WCS information to text. */
status=wcshdo(WCSHDO_safe, wcs, &nkeyrec, &wcsheader);
if(status)
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 4adda93..70ba6bc 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -51,6 +51,9 @@ double *
gal_wcs_array_from_wcsprm(struct wcsprm *wcs);
void
+gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs);
+
+void
gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy, double *radec,
size_t number, size_t width);
diff --git a/lib/wcs.c b/lib/wcs.c
index 2d053d4..4224a70 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -78,6 +78,61 @@ gal_wcs_array_from_wcsprm(struct wcsprm *wcs)
+/* According to the FITS standard, in the `PCi_j' WCS formalism, the matrix
+ elements m_{ij} are encoded in the `PCi_j' keywords and the scale
+ factors are encoded in the `CDELTi' keywords. There is also another
+ formalism (the `CDi_j' formalism) which merges the two into one
+ matrix.
+
+ However, WCSLIB's internal operations are apparently done in the `PCi_j'
+ formalism. So its outputs are also all in that format by default. When
+ the input is a `CDi_j', WCSLIB will still read the image into the
+ `PCi_j' formalism and the `CDELTi's are set to 1. This function will
+ decompose the two matrices to give a reasonable `CDELTi' and `PCi_j' in
+ such cases. */
+void
+gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
+{
+ double *ps;
+ size_t i, j;
+ int sumcond=0;
+
+ /* The correction is only needed when the matrix is internally stored
+ as PCi_j. */
+ if(wcs->altlin |= 1)
+ {
+ /* Check if the `PCi_j' and `CDELTi's have already been
+ decomposed. If all the `CDELTi's are 1.0, then most probably they
+ haven't. */
+ for(i=0; i<wcs->naxis; ++i)
+ sumcond += wcs->cdelt[i]==1.0f;
+
+ /* If all `CDELTi's were 1, then `sumcond' is equal to the number of
+ WCS axises. */
+ if(sumcond==wcs->naxis)
+ {
+ /* Get the pixel scale. */
+ ps=gal_wcs_pixel_scale_deg(wcs);
+
+ /* Correct the CDELTs. */
+ for(i=0; i<wcs->naxis; ++i)
+ wcs->cdelt[i] *= ps[i];
+
+ /* Correct the PCi_js */
+ for(i=0;i<wcs->naxis;++i)
+ for(j=0;j<wcs->naxis;++j)
+ wcs->pc[i*wcs->naxis+j] /= ps[i];
+
+ /* Clean up. */
+ free(ps);
+ }
+ }
+}
+
+
+
+
+
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 9321009: Decomposing PCi_j and CDELTi matices in output WCS,
Mohammad Akhlaghi <=