[gnuastro-commits] master 2ff65bef: Library (wcs.h): coverage func works with image that crosses RA=0

[Mohammad Akhlaghi]

branch: master
commit 2ff65befb9e07daca58a31d1b75eecd5dc76f1d8
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (wcs.h): coverage func works with image that crosses RA=0
    
    Until now, when the input image crossed the RA=0 hour circle, the outputs
    of the 'gal_wcs_coverage' function were unreasonable (for example the width
    for a ~1 degree image would be about 359 degres!).
    
    With this commit, the function finds the dimension that may have this
    problem (usually the RA dimension, but can also be the galactic or ecliptic
    longitudes for example), checks if the width is larger than 180 degrees,
    and if so, makes the necessary corrections to report the proper width and
    minimum and maximum coordinates.
    
    In the process, I also noticed an "may be used uninitialize" compiler
    warning for the function pointer 'warp_pixel_perimeter' within the
    'warp_pixelarea_onthread' function of 'lib/warp.c'. So it has been
    initialized to NULL to have a clean build with no un-necessary warnings.
    
    This bug was reported by Irene Pintos Castro.
    
    This fixes bug #63257.
---
 NEWS                         |  3 ++
 doc/announce-acknowledge.txt |  1 +
 doc/gnuastro.texi            |  9 ++++++
 lib/warp.c                   |  4 +--
 lib/wcs.c                    | 72 ++++++++++++++++++++++++++++++++++++++------
 5 files changed, 78 insertions(+), 11 deletions(-)

diff --git a/NEWS b/NEWS
index 009a9d86..ed51e328 100644
--- a/NEWS
+++ b/NEWS
@@ -315,6 +315,9 @@ See the end of the file for license conditions.
               being NaN, reported by Nafise Sedighi.
   bug #63207: Match crashes when one input has no rows. Found by Sepideh
               Eskandarlou.
+  bug #63257: Fits program's '--skycoverage' gives unreasonable outputs
+              when image crosses the RA=0 hour circle. Found by Irene
+              Pintos Castro.
 
 
 
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 53b8b205..ff0d1ebb 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -6,6 +6,7 @@ Sepideh Eskandarlou
 Giulia Golini
 Raul Infante-Sainz
 Teet Kuutma
+Irene Pintos Castro
 Nafise Sedighi
 Richard Stallman
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 019cd1ca..1bd2c8ad 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -11975,6 +11975,10 @@ $ astfits castor.fits --skycoverage --quiet
 Note that this is a simple rectangle (cube in 3D) definition, so if the image 
is rotated in relation to the celestial coordinates a general polygon is 
necessary to exactly describe the coverage.
 Hence when there is rotation, the reported area will be larger than the actual 
area containing data, you can visually see the area with the 
@option{--pixelareaonwcs} option of @ref{Fits}.
 
+Currently this option only supports images that are less than 180 degrees in 
width (which is usually the case!).
+This requirement has been necessary to account for images that cross the RA=0 
hour circle on the sky.
+Please get in touch with us at @url{mailto:bug-gnuastro@@gnu.org} if you have 
an image that is larger than 180 degrees so we try to find a solution based on 
need.
+
 @item --datasum
 @cindex @code{DATASUM}: FITS keyword
 Calculate and print the given HDU's "datasum" to stdout.
@@ -34211,6 +34215,11 @@ In other cases, this function will return a NaN.
 @deftypefun int gal_wcs_coverage (char @code{*filename}, char @code{*hdu}, 
size_t @code{*ondim}, double @code{**ocenter}, double @code{**owidth}, double 
@code{**omin}, double @code{**omax})
 Find the sky coverage of the image HDU (@code{hdu}) within @file{filename}.
 The number of dimensions is written into @code{ndim}, and space for the 
various output arrays is internally allocated and filled with the respective 
values.
+Therefore you need to free them afterwards.
+
+Currently this function only supports images that are less than 180 degrees in 
width (which is usually the case!).
+This requirement has been necessary to account for images that cross the RA=0 
hour circle on the sky.
+Please get in touch with us at @url{mailto:bug-gnuastro@@gnu.org} if you have 
an image that is larger than 180 degrees so we try to find a solution based on 
need.
 @end deftypefun
 
 @deftypefun {gal_data_t *} gal_wcs_world_to_img (gal_data_t @code{*coords}, 
struct wcsprm @code{*wcs}, int @code{inplace})
diff --git a/lib/warp.c b/lib/warp.c
index adb23af9..9a9ea439 100644
--- a/lib/warp.c
+++ b/lib/warp.c
@@ -1115,8 +1115,8 @@ warp_pixelarea_onthread(void *inparam)
 
   /* Low-level variables. */
   size_t i, ind;
-  double area, *ocrn, *outputarr=wa->output->array;
-  double *(*warp_pixel_perimeter)(gal_warp_wcsalign_t *, size_t);
+  double area, *ocrn=NULL, *outputarr=wa->output->array;
+  double *(*warp_pixel_perimeter)(gal_warp_wcsalign_t *, size_t)=NULL;
 
   /* Call the correct function based on the output image orientation. */
   if( wa->isccw==1 )
diff --git a/lib/wcs.c b/lib/wcs.c
index 1f98553a..9332e6c2 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -2040,9 +2040,10 @@ gal_wcs_coverage(char *filename, char *hdu, size_t 
*ondim,
   struct wcsprm *wcs;
   int nwcs=0, type, status=0;
   char *name=NULL, *unit=NULL;
-  gal_data_t *tmp, *coords=NULL;
+  gal_data_t *fc, *tmp, *coords=NULL;
   size_t i, ndim, *dsize=NULL, numrows;
-  double *x=NULL, *y=NULL, *z=NULL, *min, *max, *center, *width;
+  size_t fullcircledim=GAL_BLANK_SIZE_T;
+  double *x=NULL, *y=NULL, *z=NULL, *min, *max, *fcarr, *center, *width;
 
   /* Read the desired WCS (note that the linear matrix is irrelevant here,
      we'll just select PC because its the default WCS mode. */
@@ -2161,20 +2162,73 @@ gal_wcs_coverage(char *filename, char *hdu, size_t 
*ondim,
     }
 
   /* Write the center and width. */
+  width[0]=max[0]-min[0];
+  width[1]=max[1]-min[1];
   switch(ndim)
     {
     case 2:
-      center[0]=x[4];         center[1]=y[4];
-      width[0]=max[0]-min[0]; width[1]=max[1]-min[1];
+      center[0]=x[4];  center[1]=y[4];
       break;
     case 3:
-      center[0]=x[8];         center[1]=y[8];         center[2]=z[8];
-      width[0]=max[0]-min[0]; width[1]=max[1]-min[1]; width[2]=max[2]-min[2];
+      width[2]=max[2]-min[2];
+      center[0]=x[8]; center[1]=y[8]; center[2]=z[8];
       break;
     default:
-      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to solve the "
-            "problem. The value %zu is not a recognized dimension", __func__,
-            PACKAGE_BUGREPORT, ndim);
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to solve "
+            "the problem. The value %zu is not a recognized dimension",
+            __func__, PACKAGE_BUGREPORT, ndim);
+    }
+
+  /* It may happen that the image passes the RA=0 hour circle. In that
+     case, the width along the first WCS dimension (RA) will be larger than
+     180. Therefore, we need to re-calculate the minimum and maximums. An
+     example dataset is available in https://savannah.gnu.org/bugs/?63257
+     (Gnuastro bug #63257):
+
+        $ jplusdr2url=http://archive.cefca.es/catalogues/vo/siap/jplus-dr2
+        $ wget $jplusdr2url/get_fits?id=71811 -Ojplus-dr2-71811.fits.fz
+
+     But first, we need to find the dimension that has a full-circle
+     dimension: this is usually the first WCS dimension, but it may happen
+     that it isn't. Also, note that the WCS may not be celestial/spherical
+     at all (hence why we are only doing this for pre-defined celestial
+     'CTYPE's).*/
+  for(i=0;i<ndim;++i)
+    if(   strncmp(wcs->ctype[i], "RA-",   3)==0 /* Equatorial */
+       || strncmp(wcs->ctype[i], "GLON-", 5)==0 /* Galactic */
+       || strncmp(wcs->ctype[i], "SLON-", 5)==0 /* Super Galactic */
+       || strncmp(wcs->ctype[i], "ELON-", 5)==0 /* Ecliptic */ )
+      { fullcircledim=i; break; }
+  if( fullcircledim!=GAL_BLANK_SIZE_T && width[fullcircledim]>180.0f )
+    {
+      /* Set the proper pointer to the coordinate that should be
+         checked. */
+      switch(fullcircledim)
+        {
+        case 0: fc=coords;             break;
+        case 1: fc=coords->next;       break;
+        case 2: fc=coords->next->next; break;
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+                "'%s' to fix the problem. The value '%zu' isn't "
+                "recognized for 'fullcircledim'", __func__,
+                PACKAGE_BUGREPORT, fullcircledim);
+        }
+
+      /* Shift all the coordinates after 0, by 360, while keeping those
+         that are close to 360 untouched, then find the minimum and
+         maximums. */
+      fcarr=fc->array;
+      for(i=0;i<coords->size;++i) if(fcarr[i]<180) fcarr[i]+=360.0f;
+      tmp=gal_statistics_minimum(fc);
+      min[fullcircledim] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+      tmp=gal_statistics_maximum(fc);
+      max[fullcircledim] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+
+      /* Calculate the correct width, maximum is guarateed to be larger
+         than 360, so subtract 360 from it. */
+      width[fullcircledim]=max[fullcircledim]-min[fullcircledim];
+      max[fullcircledim]-=360.0f;
     }
 
   /* Clean up and return success. */