[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 4b747b0: Check on floating point errors in pix
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 4b747b0: Check on floating point errors in pixel scale calculation |
Date: |
Thu, 21 Sep 2017 20:32:17 -0400 (EDT) |
branch: master
commit 4b747b037b280b1c5737bff227a3bb9a04375d21
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Check on floating point errors in pixel scale calculation
See the added error description for the full discussion: in short, some
programs include floating point errors in their WCS FITS header keywords
and these can bias the pixel scale measurements (and produce an error
message sometimes). So with this commit, a warning will now be printed to
guide the users on the (possible) problem and how to solve it.
This was brought up in a discussion with Leindert Boogaard.
---
THANKS | 1 +
doc/announce-acknowledge.txt | 1 +
lib/wcs.c | 59 ++++++++++++++++++++++++++++++++++++++++++--
3 files changed, 59 insertions(+), 2 deletions(-)
diff --git a/THANKS b/THANKS
index 5f34189..1afb512 100644
--- a/THANKS
+++ b/THANKS
@@ -18,6 +18,7 @@ support in Gnuastro. The list is ordered alphabetically (by
family name).
Marjan Akbari address@hidden
Karl Berry address@hidden
Roland Bacon address@hidden
+ Leindert Boogaard address@hidden
Nicolas Bouché address@hidden
Fernando Buitrago address@hidden
Adrian Bunk address@hidden
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 880c69c..4e2bfea 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,5 +1,6 @@
This file is meant to keep the names of the people who's help must be
acknowledged in the next release.
+Leindert Boogaard
Takashi Ichikawa
David Valls-Gabaud
diff --git a/lib/wcs.c b/lib/wcs.c
index 1fe18da..5b438eb 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -438,10 +438,10 @@ double *
gal_wcs_pixel_scale(struct wcsprm *wcs)
{
gsl_vector S;
- double *a, *out;
- int permute_set;
gsl_matrix A, V;
size_t i, j, n=wcs->naxis;
+ double *a, *out, maxrow, minrow;
+ int permute_set, warning_printed;
double *v=gal_data_malloc_array(GAL_TYPE_FLOAT64, n*n, __func__, "v");
size_t *permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, n, __func__,
"permutation");
@@ -455,6 +455,61 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
a=gal_wcs_warp_matrix(wcs);
+ /* To avoid confusing issues with floating point errors being written in
+ the non-diagonal elements of the FITS header PC or CD matrices, we
+ need to check if the minimum and maximum values in each row are not
+ several orders of magnitude apart.
+
+ Note that in some cases (for example a spectrum), one axis might be in
+ degrees (value around 1e-5) and the other in angestroms (value around
+ 1e-10). So we can't look at the minimum and maximum of the whole
+ matrix. However, in such cases, people will probably not warp/rotate
+ the image to mix the coordinates. So the important thing to check is
+ the minimum and maximum (non-zero) values in each row. */
+ warning_printed=0;
+ for(i=0;i<n;++i)
+ {
+ /* Find the minimum and maximum values in each row. */
+ minrow=FLT_MAX;
+ maxrow=-FLT_MAX;
+ for(j=0;j<n;++j)
+ if(a[i*n+j]!=0.0) /* We aren't concerned with 0 valued elements. */
+ {
+ /* We have to use absolutes because in cases like RA, the
+ diagonal values in different rows may have different signs*/
+ if(fabs(a[i*n+j])<minrow) minrow=fabs(a[i*n+j]);
+ if(fabs(a[i*n+j])>maxrow) maxrow=fabs(a[i*n+j]);
+ }
+
+ /* Do the check, print warning and make correction. */
+ if(maxrow!=minrow && maxrow/minrow>1e4 && warning_printed==0)
+ {
+ fprintf(stderr, "\nWARNING: The input WCS matrix (possibly taken "
+ "from the FITS header keywords starting with `CD' or `PC') "
+ "contains values with very different scales (more than "
+ "10^4 different). This is probably due to floating point "
+ "errors. These values might bias the pixel scale (and "
+ "subsequent) calculations.\n\n"
+ "You can see the respective matrix with one of the "
+ "following two commands (depending on how the FITS file "
+ "was written). Recall that if the desired extension/HDU "
+ "isn't the default, you can choose it with the `--hdu' "
+ "(or `-h') option before the `|' sign in these commands."
+ "\n\n"
+ " $ astfits file.fits -p | grep 'PC._.'\n"
+ " $ astfits file.fits -p | grep 'CD._.'\n\n"
+ "You can delete the ones with obvious floating point "
+ "error values using the following command (assuming you "
+ "want to delete `CD1_2' and `CD2_1') and rerun the your "
+ "command to remove this warning message and possibly "
+ "correct errors that might be reported due to it.\n\n"
+ " $ astfits file.fits --delete=CD1_2 --delete=CD2_1\n\n"
+ );
+ warning_printed=1;
+ }
+ }
+
+
/* Fill in the necessary GSL vector and Matrix structures. */
S.size=n; S.stride=1; S.data=pixscale->array;
V.size1=n; V.size2=n; V.tda=n; V.data=v;
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 4b747b0: Check on floating point errors in pixel scale calculation,
Mohammad Akhlaghi <=