[gnuastro-commits] master fac42e64: Library (wcs.h): conversion functions always return, NaN when bad

[Mohammad Akhlaghi]

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

    Library (wcs.h): conversion functions always return, NaN when bad
    
    Until now, if only a single one of the input coordiantes was not valid (or
    had any other kind of WCSLIB related error), the two 'gal_wcs_img_to_word'
    and 'gal_wcs_world_to_img' functions would abort the program. However, it
    often happens that the user doesn't care about the bad points and can later
    remove them!
    
    With this commit, we do a coordinate-by-coordinate check and set the bad
    coordinates to NaN. This lets the users select if they want to keep or
    remove that particular coordinate.
---
 doc/gnuastro.texi | 42 ++++++++++++++++--------------------------
 lib/wcs.c         | 41 ++++++++++++++++++++++-------------------
 2 files changed, 38 insertions(+), 45 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 388c42e6..99882287 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -32957,19 +32957,16 @@ points is calculated with the equation below.
 
 @dispmath {\cos(d)=\sin(d_1)\sin(d_2)+\cos(d_1)\cos(d_2)\cos(r_1-r_2)}
 
-However, since the pixel scales are usually very small numbers, this
-function will not use that direct formula. It will be use the
-@url{https://en.wikipedia.org/wiki/Haversine_formula, Haversine formula}
-which is better considering floating point errors:
+However, since the pixel scales are usually very small numbers, this function 
will not use that direct formula.
+It will be use the @url{https://en.wikipedia.org/wiki/Haversine_formula, 
Haversine formula} which is better considering floating point errors:
 
 @dispmath{{\sin^2(d)\over 2}=\sin^2\left( {d_1-d_2\over 2} 
\right)+\cos(d_1)\cos(d_2)\sin^2\left( {r_1-r_2\over 2} \right)}
 @end deftypefun
 
 @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 allocated array of double precision floating point type with
-one element for each dimension. If it is not successful, this function will
-return @code{NULL}.
+Return the pixel scale for each dimension of @code{wcs} in degrees.
+The output is an allocated array of double precision floating point type with 
one element for each dimension.
+If it is not successful, this function will return @code{NULL}.
 @end deftypefun
 
 @deftypefun double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
@@ -32984,27 +32981,20 @@ The number of dimensions is written into @code{ndim}, 
and space for the various
 @end deftypefun
 
 @deftypefun {gal_data_t *} gal_wcs_world_to_img (gal_data_t @code{*coords}, 
struct wcsprm @code{*wcs}, int @code{inplace})
-Convert the linked list of world coordinates in @code{coords} to a linked
-list of image coordinates given the input WCS structure. @code{coords} must
-be a linked list of data structures of float64 (`double') type,
-see@ref{Linked lists} and @ref{List of gal_data_t}. The top (first
-popped/read) node of the linked list must be the first WCS coordinate (RA
-in an image usually) etc. Similarly, the top node of the output will be
-the first image coordinate (in the FITS standard).
-
-If @code{inplace} is zero, then the output will be a newly allocated list
-and the input list will be untouched. However, if @code{inplace} is
-non-zero, the output values will be written into the input's already
-allocated array and the returned pointer will be the same pointer to
-@code{coords} (in other words, you can ignore the returned value). Note
-that in the latter case, only the values will be changed, things like units
-or name (if present) will be untouched.
+Convert the linked list of world coordinates in @code{coords} to a linked list 
of image coordinates given the input WCS structure.
+@code{coords} must be a linked list of data structures of float64 (`double') 
type, see@ref{Linked lists} and @ref{List of gal_data_t}.
+The top (first popped/read) node of the linked list must be the first WCS 
coordinate (RA in an image usually) etc.
+Similarly, the top node of the output will be the first image coordinate (in 
the FITS standard).
+In case WCSLIB fails to convert any of the coordinates (for example the RA of 
one coordinate is given as 400!), the respective element in the output will be 
written as NaN.
+
+If @code{inplace} is zero, then the output will be a newly allocated list and 
the input list will be untouched.
+However, if @code{inplace} is non-zero, the output values will be written into 
the input's already allocated array and the returned pointer will be the same 
pointer to @code{coords} (in other words, you can ignore the returned value).
+Note that in the latter case, only the values will be changed, things like 
units or name (if present) will be untouched.
 @end deftypefun
 
 @deftypefun {gal_data_t *} gal_wcs_img_to_world (gal_data_t @code{*coords}, 
struct wcsprm @code{*wcs}, int @code{inplace})
-Convert the linked list of image coordinates in @code{coords} to a linked
-list of world coordinates given the input WCS structure. See the
-description of @code{gal_wcs_world_to_img} for more details.
+Convert the linked list of image coordinates in @code{coords} to a linked list 
of world coordinates given the input WCS structure.
+See the description of @code{gal_wcs_world_to_img} for more details.
 @end deftypefun
 
 
diff --git a/lib/wcs.c b/lib/wcs.c
index 23e489de..47be5227 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -2230,7 +2230,7 @@ wcs_convert_sanity_check_alloc(gal_data_t *coords, struct 
wcsprm *wcs,
   /* Allocate all the necessary arrays. */
   *phi    = gal_pointer_allocate( GAL_TYPE_FLOAT64, size,      0, __func__,
                                   "phi");
-  *stat   = gal_pointer_allocate( GAL_TYPE_INT32,   size,      1, __func__,
+  *stat   = gal_pointer_allocate( GAL_TYPE_INT,     size,      1, __func__,
                                   "stat");
   *theta  = gal_pointer_allocate( GAL_TYPE_FLOAT64, size,      0, __func__,
                                   "theta");
@@ -2284,14 +2284,15 @@ wcs_convert_prepare_out(gal_data_t *coords, struct 
wcsprm *wcs, int inplace)
 {
   size_t i;
   gal_data_t *out=NULL;
+
   if(inplace)
     out=coords;
   else
     for(i=0;i<wcs->naxis;++i)
       gal_list_data_add_alloc(&out, NULL, GAL_TYPE_FLOAT64, 1,
                               &coords->size, NULL, 0, coords->minmapsize,
-                              coords->quietmmap, wcs->ctype[i], wcs->cunit[i],
-                              NULL);
+                              coords->quietmmap, wcs->ctype[i],
+                              wcs->cunit[i], NULL);
   return out;
 }
 
@@ -2308,7 +2309,7 @@ gal_data_t *
 gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm *wcs, int inplace)
 {
   gal_data_t *out;
-  int status, *stat=NULL, ncoord=coords->size, nelem;
+  int *stat=NULL, ncoord=coords->size, nelem;
   double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
 
   /* It can happen that the input datasets are empty. In this case, simply
@@ -2321,6 +2322,7 @@ gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm 
*wcs, int inplace)
                  "'inplace' is not called", __func__, PACKAGE_BUGREPORT);
     }
 
+
   /* Some sanity checks. */
   wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi,
                                  &theta, &world, &pixcrd, &imgcrd);
@@ -2332,12 +2334,10 @@ gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm 
*wcs, int inplace)
   wcs_convert_list_to_from_array(coords, world, stat, wcs->naxis, 0);
 
 
-  /* Use WCSLIB's wcss2p for the conversion. */
-  status=wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd,
-                pixcrd, stat);
-  if(status)
-    error(EXIT_FAILURE, 0, "%s: wcss2p ERROR %d: %s", __func__, status,
-          wcs_errmsg[status]);
+  /* Use WCSLIB's wcss2p for the conversion. We are ignoring the over-all
+     status here, because later we will use the 'stat' array to set all bad
+     coordinates to NaN. */
+  wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd, pixcrd, stat);
 
 
   /* For a sanity check.
@@ -2345,9 +2345,13 @@ gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm 
*wcs, int inplace)
     size_t i;
     printf("\n\n%s sanity check:\n", __func__);
     for(i=0;i<coords->size;++i)
-      printf("(%g, %g, %g) --> (%g, %g, %g), [stat: %d]\n",
-              world[i*3],  world[i*3+1 ], world[i*3+2],
-             pixcrd[i*3], pixcrd[i*3+1], pixcrd[i*3+2], stat[i]);
+      printf("(%g, %g) --> phi:%g, theta: %g --> (%g, %g) --> (%g, %g) "
+             "[stat: %d]\n",
+             world[i*2],  world[i*2+1],
+             phi[i*2],    theta[i*2],
+             imgcrd[i*2], imgcrd[i*2+1],
+             pixcrd[i*2], pixcrd[i*2+1], stat[i]);
+    printf("stat: %d\n", stat[4]);
   }
   */
 
@@ -2382,7 +2386,7 @@ gal_data_t *
 gal_wcs_img_to_world(gal_data_t *coords, struct wcsprm *wcs, int inplace)
 {
   gal_data_t *out;
-  int status, *stat=NULL, ncoord=coords->size, nelem;
+  int *stat=NULL, ncoord=coords->size, nelem;
   double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
 
   /* Some sanity checks. */
@@ -2396,11 +2400,10 @@ gal_wcs_img_to_world(gal_data_t *coords, struct wcsprm 
*wcs, int inplace)
   wcs_convert_list_to_from_array(coords, pixcrd, stat, wcs->naxis, 0);
 
 
-  /* Use WCSLIB's wcsp2s for the conversion. */
-  status=wcsp2s(wcs, ncoord, nelem, pixcrd, imgcrd, phi, theta, world, stat);
-  if(status)
-    error(EXIT_FAILURE, 0, "%s: wcsp2s ERROR %d: %s", __func__, status,
-          wcs_errmsg[status]);
+  /* Use WCSLIB's wcsp2s for the conversion. We are ignoring the over-all
+     status here, because later we will use the 'stat' array to set all bad
+     coordinates to NaN.*/
+  wcsp2s(wcs, ncoord, nelem, pixcrd, imgcrd, phi, theta, world, stat);
 
 
   /* For a check.