[gnuastro-commits] master 4456ad3: Library (WCS): will not crash with singular WCS matrix

[Mohammad Akhlaghi]

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

    Library (WCS): will not crash with singular WCS matrix
    
    Until now, there was no check for a singular WCS matrix (where the
    determinant is zero). But in one dataset, we noticed that this situation
    has occurrred, so there were bad warnings and crashs when running
    NoiseChisel.
    
    With this commit, in such situations, the pixel scale (which caused those
    errors) will not be calculated and the higher-level functions that used it
    are also able to avoid the problem. Ofcourse a problem may arise in the
    interpretation of the WCS and in that case, it should indeed crash. But in
    writing a dataset, it shouldn't.
    
    Also, while testing with the dataset (in the middle of the checks above), I
    noticed that there are no checks to avoid a crash in Fits when called with
    `--copykeys', but without any output filename or HDU.
---
 bin/crop/wcsmode.c |  3 +++
 bin/fits/args.h    |  4 +++-
 bin/fits/ui.c      | 17 +++++++++++------
 bin/warp/ui.c      |  3 +++
 doc/gnuastro.texi  |  5 +++--
 lib/wcs.c          | 30 +++++++++++++++++++++++-------
 6 files changed, 46 insertions(+), 16 deletions(-)

diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index be5101e..04ea1d0 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -101,6 +101,9 @@ wcsmode_check_prepare(struct cropparams *p, struct 
inputimgs *img)
      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(pixscale==NULL)
+    error(EXIT_FAILURE, 0, "the pixel scale couldn't be deduced from the "
+          "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 "
diff --git a/bin/fits/args.h b/bin/fits/args.h
index 2a1360b..5f3d899 100644
--- a/bin/fits/args.h
+++ b/bin/fits/args.h
@@ -92,6 +92,8 @@ struct argp_option program_options[] =
 
 
 
+
+
     {
       0, 0, 0, 0,
       "Keywords (in one HDU):",
@@ -230,7 +232,7 @@ struct argp_option program_options[] =
     {
       "copykeys",
       UI_KEY_COPYKEYS,
-      "STR",
+      "INT:INT",
       0,
       "Range of keywords to copy to output HDU.",
       UI_GROUP_KEYWORD,
diff --git a/bin/fits/ui.c b/bin/fits/ui.c
index c0c14d7..1d69fd4 100644
--- a/bin/fits/ui.c
+++ b/bin/fits/ui.c
@@ -221,9 +221,15 @@ ui_check_copykeys(struct fitsparams *p)
 {
   long read;
   char *tailptr;
-  size_t group=0;
+  /* size_t group=0; */
   char forl='f', *pt=p->copykeys;
 
+  /* For copykeys, an output filename is mandatory. */
+  if(p->cp.output==NULL || p->outhdu==NULL)
+    error(EXIT_FAILURE, 0, "an output FITS extension (in an existing "
+          "FITS file, specified with the `--output' and `--outhdu') are "
+          "mandatory for running `--copykeys'");
+
   /* Initialize the values. */
   p->copykeysrange[0]=p->copykeysrange[1]=GAL_BLANK_LONG;
 
@@ -251,20 +257,19 @@ ui_check_copykeys(struct fitsparams *p)
         case '6': case '7': case '8': case '9': case '-':
           break;
 
-        /* An un-recognized character should crash the program. */
+        /* Identifier for next group of ranges. However, For the time
+           being, we just support one group. So we are commenting the break
+           here for it to follow onto default.
         case ',':
           ++group;
           forl='f';
           ++pt;
-
-          /* For the time being, we just support one group. So we are
-             commenting the break here for it to follow onto default.
           break;
           */
         default:
           error(EXIT_FAILURE, 0, "value to `--copykeys' must only contain "
                 "integer numbers and these special characters between them: "
-                "`,', `:', `*' when necessary. But it is `%s' (the first "
+                "`:' when necessary. But it is `%s' (the first "
                 "non-acceptable character is `%c').\n", p->copykeys, *pt);
           break;
         }
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 3e2d8e7..3c09b7d 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -348,6 +348,9 @@ ui_check_options_and_arguments(struct warpparams *p)
       if(p->input->wcs)
         {
           p->pixelscale=gal_wcs_pixel_scale(p->input->wcs);
+          if(p->pixelscale==NULL)
+            error(EXIT_FAILURE, 0, "%s (hdu %s): the pixel scale couldn't "
+                  "be deduced from the WCS.", p->inputname, p->cp.hdu);
           p->inwcsmatrix=gal_wcs_warp_matrix(p->input->wcs);
         }
     }
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c05151b..684636f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -27347,8 +27347,9 @@ which is better considering floating point errors:
 
 @deftypefun {double *} gal_wcs_pixel_scale (struct wcsprm @code{*wcs})
 Return the pixel scale for each dimension of @code{wcs} in degrees. The
-output is an array of double precision floating point type with one element
-for each dimension.
+output is an allocated array of double precision floating point type with
+one element for each dimension. If its not successful, this function will
+return @code{NULL}.
 @end deftypefun
 
 @deftypefun double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
diff --git a/lib/wcs.c b/lib/wcs.c
index 23b5389..9e3c00c 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -564,11 +564,15 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
   double *ps;
   size_t i, j;
 
+  /* If there is on WCS, then don't do anything. */
+  if(wcs==NULL) return;
+
   /* The correction is only needed when the PC matrix is filled. */
   if(wcs->pc)
     {
       /* Get the pixel scale. */
       ps=gal_wcs_pixel_scale(wcs);
+      if(ps==NULL) return;
 
       /* The PC matrix and the CDELT elements might both contain scale
          elements (during processing the scalings might be added only to PC
@@ -636,6 +640,7 @@ gal_wcs_angular_distance_deg(double r1, double d1, double 
r2, double d2)
 
 
 
+
 /* Return the pixel scale of the dataset in units of the WCS. */
 double *
 gal_wcs_pixel_scale(struct wcsprm *wcs)
@@ -643,20 +648,30 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
   gsl_vector S;
   gsl_matrix A, V;
   int warning_printed;
-  size_t i, j, maxj, n=wcs->naxis;
-  double jvmax, *a, *out, maxrow, minrow;
-  double *v=gal_pointer_allocate(GAL_TYPE_FLOAT64, n*n, 0, __func__, "v");
-  size_t *permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, n, 0, __func__,
-                                           "permutation");
-  gal_data_t *pixscale=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &n, NULL,
-                                      0, -1, NULL, NULL, NULL);
+  gal_data_t *pixscale;
+  size_t i, j, n, maxj, *permutation;
+  double jvmax, *a, *out, *v, maxrow, minrow;
 
+  /* Only continue if a WCS exists. */
+  if(wcs==NULL) return NULL;
 
   /* Write the full WCS rotation matrix into an array, irrespective of what
      style it was stored in the wcsprm structure (`PCi_j' style or `CDi_j'
      style). */
   a=gal_wcs_warp_matrix(wcs);
 
+  /* A small sanity check (this won't work on a singular matrix, can happen
+     in FITS WCSs!). In this case, we should return NULL.*/
+  n=wcs->naxis;
+  for(i=0;i<n;++i) {if(a[i*n+i]==0.0f) return NULL;}
+
+  /* Now that everything is good, we can allocate the necessary memory. */
+  v=gal_pointer_allocate(GAL_TYPE_FLOAT64, n*n, 0, __func__, "v");
+  permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, n, 0, __func__,
+                                   "permutation");
+  pixscale=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &n, NULL,
+                          0, -1, NULL, NULL, NULL);
+
 
   /* To avoid confusing issues with floating point errors being written in
      the non-diagonal elements of the FITS header PC or CD matrices, we
@@ -792,6 +807,7 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
 
   /* Get the pixel scales along each axis in degrees, then multiply. */
   pixscale=gal_wcs_pixel_scale(wcs);
+  if(pixscale==NULL) return NAN;
 
   /* Clean up and return the result. */
   out = pixscale[0] * pixscale[1] * 3600.0f * 3600.0f;