[gnuastro-commits] master 9321009: Decomposing PCi_j and CDELTi matices in output WCS

[Mohammad Akhlaghi]

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);
+        }
+    }
+}
+
+
+
+
+