[gnuastro-commits] master e536c79: Always decompose PC and CDELT elements

[Mohammad Akhlaghi]

branch: master
commit e536c797a874f56fd15e6ecac9a69b44f3b919fd
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Always decompose PC and CDELT elements
    
    Previously the decomposition of the `PCi_j' matrix and the `CDELTi' vector
    would only be done when the elements of the latter were only 1. Either from
    the input, or during the processing, the final values of `PCi_j' might
    contain scale factors in them, independent of the `CDELTi' vector. For
    example this is what ImageWarp currently does (it just modifies the `PCi_j'
    matrix).
    
    This was visible in the headers of `convolve_spatial_warped.fits' (after
    running `make check'): The `PCi_j' diagonal matrix were a multiple of 5 and
    the `CDELTi' was the original (input resolution).
    
    So with this commit, the decomposition of `PCi_j' and `CDELTi' is done all
    cases when a `PCi_j' matrix is present: the pixel scale is found, they are
    unified, then the pixel scale in each dimension is put into `CDELTi' and
    divided by `PCi_j'.
---
 lib/wcs.c |   51 ++++++++++++++++++++++++++++-----------------------
 1 file changed, 28 insertions(+), 23 deletions(-)

diff --git a/lib/wcs.c b/lib/wcs.c
index 4224a70..a23c6c0 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -95,37 +95,42 @@ 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);
+      /* 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];
+      /* The PC matrix and the CDELT elements might both contain scale
+         elements (during processing the scalings might be added only to PC
+         elements for example). So to be safe, we first multiply them into
+         one matrix. */
+      for(i=0;i<wcs->naxis;++i)
+        for(j=0;j<wcs->naxis;++j)
+          wcs->pc[i*wcs->naxis+j] *= wcs->cdelt[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];
+      /* Set the CDELTs. */
+      for(i=0; i<wcs->naxis; ++i)
+        wcs->cdelt[i] = ps[i];
 
-          /* Clean up. */
-          free(ps);
-        }
+      /* 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);
+
+      /* According to the `wcslib/wcs.h' header: "In particular, wcsset()
+         resets wcsprm::cdelt to unity if CDi_ja is present (and no
+         PCi_ja).". So apparently, when the input is a `CDi_j', it might
+         expect the `CDELTi' elements to be 1.0. But we have changed that
+         here, so we will correct the `altlin' element of the WCS structure
+         to make sure that WCSLIB only looks into the `PCi_j' and `CDELTi'
+         and makes no assumptioins about `CDELTi'. */
+      wcs->altlin=1;
     }
 }