[gnuastro-commits] master 0dec9ac 110/125: Allow for floating point errors in Crop comparisons

[Mohammad Akhlaghi]

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

    Allow for floating point errors in Crop comparisons
    
    When the input image is not aligned with the WCS or the pixel scales along
    each dimension are different, in WCS mode, Crop will complain and abort
    (because it currently uses these assumptions). Until now this check was
    done exactly, so even a small difference of, for example 1e-10 (which is
    clearly due to floating point errors) would cause an error. So with this
    commit, very small relative differences (less than 1e-6) are acceptable and
    Crop won't complain.
    
    In this commit, short explanation that was also added to the `--polygon'
    option of Crop so users can use SAO ds9 (a FITS viewer) for defining the
    polygon instead of having to type the values manually.
    
    Both these suggestions were kindly made by Guillaume Mahler, so his name
    has been added to the `THANKS' file.
---
 THANKS             |  1 +
 bin/crop/crop.c    |  4 ++++
 bin/crop/wcsmode.c | 33 ++++++++++++++++++---------------
 doc/gnuastro.texi  | 43 +++++++++++++++++++++++++++++++++++++------
 4 files changed, 60 insertions(+), 21 deletions(-)

diff --git a/THANKS b/THANKS
index b1d13c7..a28ae7c 100644
--- a/THANKS
+++ b/THANKS
@@ -23,6 +23,7 @@ support in Gnuastro. The list is ordered alphabetically.
     Brandon Invergo                      address@hidden
     Lee Kelvin                           address@hidden
     Mohammad-Reza Khellat                address@hidden
+    Guillaume Mahler                     address@hidden
     Alan Lefor                           address@hidden
     Francesco Montanari                  address@hidden
     William Pence                        address@hidden
diff --git a/bin/crop/crop.c b/bin/crop/crop.c
index 12a708e..2de990d 100644
--- a/bin/crop/crop.c
+++ b/bin/crop/crop.c
@@ -308,6 +308,10 @@ wcsmodecrop(void *inparam)
           }
       while ( ++(crp->in_ind) < p->numin );
 
+      /* Correct in_ind. The loop above went until `in_ind' is one more
+         than the index for the last input image (that is how it exited the
+         loop). But `crp->in_ind' is needed later, so correct it here. */
+      --crp->in_ind;
 
       /* Check the final output: */
       if(crp->numimg)
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index 3173ccf..77f3684 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -56,15 +56,18 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs 
*img)
   double imgcrd[8], phi[4], theta[4], pixcrd[8];
 
 
-  /* Check if the image is aligned correctly and fits the
-     resolution of other images. */
-  if(wcs->pc[1]!=0.0f || wcs->pc[2]!=0.0f)
+  /* Check if the image is aligned with the WCS coordinates. Note that
+     because of small floating point errors, some programs might still keep
+     very small values in the off-diagonal matrix elements. */
+  if( fabs(wcs->pc[1]/wcs->pc[3])>1e-6 || fabs(wcs->pc[2]/wcs->pc[3])>1e-6 )
     error(EXIT_FAILURE, 0, "%s: HDU %s: is not aligned to the "
           "celestial coordinates. The first FITS axis should be "
           "along the Right Ascension and the second FITS axis "
-          "should be along the declination. You should rotate "
-          "(interpolate) the images with other software",
-          img->name, p->cp.hdu);
+          "should be along the declination.\n\n"
+          "Gnuastro's Warp program can align it with the following "
+          "command:\n\n"
+          "    $ astwarp %s --hdu=%s --align\n",
+          img->name, p->cp.hdu, img->name, p->cp.hdu);
   if(wcs->pc[0]>0)
     error(EXIT_FAILURE, 0, "%s: HDU %s: An increase in the first "
           "FITS axis pixel coordinates should be a decrese in the "
@@ -78,19 +81,20 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs 
*img)
           img->name, p->cp.hdu);
 
 
-  /* Since we are dealing with very accurate values, a multiplication
-     by -1 might cause a floating point error. So we have to account
-     for the floating point error. */
-  if(-1.0f*wcs->pc[0]<wcs->pc[3]-1e-15 || -1.0f*wcs->pc[0]>wcs->pc[3]+1e-15)
+  /* Check the image resolution (if the pixels are actually a square), then
+     compare the resolution with the other input images. Due to floating
+     point errors, some very 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_deg(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 "
-          "is %f arcseconds/pixel, while the second is %f",
-          img->name, p->cp.hdu, 3600*-1.0f*wcs->pc[0],
-          3600*wcs->pc[3]);
+          "is %.15g deg/pixel, while the second is %.15g",
+          img->name, p->cp.hdu, pixscale[1], pixscale[0]);
   if(p->res==0.0f)
     {
       /* Set the resolution of the image. */
-      pixscale=gal_wcs_pixel_scale_deg(wcs);
       p->res=pixscale[0];
       free(pixscale);
 
@@ -114,7 +118,6 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs 
*img)
     }
   else
     {
-      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 "
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0a8c8de..9717e7f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -7853,14 +7853,14 @@ Crop box parameters:
 
 @item -x FLT
 @itemx --xc=FLT
-X axis (first image axis, horizontal when viewed in SAO DS9) position of
+X axis (first image axis, horizontal when viewed in SAO ds9) position of
 crop box center. Along with @option{--yc}, this only produces one crop in
 each run. The width of the square crop (in pixels) is the value to
 @option{--iwidth}.
 
 @item -y FLT
 @itemx --yc=FLT
-Y axis (second image axis, vertical when viewed in SAO DS9) position of
+Y axis (second image axis, vertical when viewed in SAO ds9) position of
 crop box center. See @option{--xc}.
 
 @item -r FLT
@@ -7929,10 +7929,41 @@ have set the operating mode to WCS mode in your nearest 
configuration file
 images/tiles and Crop will stitch them to produce this cropped region:
 
 @example
-$ astcrop --mode=wcs desired-filter-image(s).fits            \
-   --polygon="53.187414,-27.779152 : 53.159507,-27.759633 :  \
+$ astcrop --mode=wcs desired-filter-image(s).fits           \
+   --polygon="53.187414,-27.779152 : 53.159507,-27.759633 : \
               53.134517,-27.787144 : 53.161906,-27.807208"
 @end example
address@hidden SAO ds9
+In other cases, you have an image and want to define the polygon yourself
+(it isn't already published like the example above). As the number of
+vertices increases, checking the vertice coordinates on a FITS viewer (for
+example SAO ds9) and typing them in one by one can be very tedious and
+prone to typo errors.
+
+You can take the following steps to avoid the frustration and possible
+typos: Open the image with ds9 and activate its ``region'' mode with
address@hidden@click{}Region}. Then define the region as a polygon
+with @address@hidden@click{}Polygon}. Click on the
+approximate center of the region you want and a small square will
+appear. By clicking on the vertices of the square you can shrink or expand
+it, clicking and dragging anywhere on the edges will enable you to define a
+new vertice. After the region has been nicely defined, save it as a file
+with @address@hidden Regions}. You can then select the
+name and address of the output file, keep the format as @command{REG} and
+press ``OK''. In the next window, keep format as ``ds9'' and ``Coordinate
+System'' as ``fk5''. A plain text file (let's call it @file{ds9.reg}) is
+now created.
+
+You can now convert this plain text file to Crop's polygon format with this
+command (when typing on the command-line, ignore the address@hidden'' at the 
end
+of the first and second lines along with the extra spaces, these are only
+for nice printing):
+
address@hidden
+$ v=$(awk 'NR==4' ds9.reg | sed -e's/polygon(//'        \
+           -e's/\([^,]*,[^,]*\),/\1:/g' -e's/)//' )
+$ astcrop --mode=wcs image.fits --polygon=$v
address@hidden example
 
 @item --outpolygon
 Keep all the regions outside the polygon and mask the inner ones with blank
@@ -10601,7 +10632,7 @@ the examples above as a demonstration.
 @cindex FITS standard
 Based on the FITS standard, integer values are assigned to the center of a
 pixel and the coordinate [1.0, 1.0] is the center of the first pixel
-(bottom left of the image when viewed in SAO DS9). So the coordinate center
+(bottom left of the image when viewed in SAO ds9). So the coordinate center
 [0.0, 0.0] is half a pixel away (in each axis) from the bottom left vertice
 of the first pixel. The resampling that is done in Warp (see
 @ref{Resampling}) is done on the coordinate axises and thus directly
@@ -10716,7 +10747,7 @@ ignored.
 @item -c
 @itemx --centeroncorer
 Put the center of coordinates on the corner of the first (bottom-left when
-viewed in SAO DS9) pixel. This option is applied after the final warping
+viewed in SAO ds9) pixel. This option is applied after the final warping
 matrix has been finalized: either through modular warpings or the raw
 matrix. See the explanation above for coordinates in the FITS standard to
 better understand this option and when it should be used.