gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 4ddc277 1/2: Correction in alignment and getti


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 4ddc277 1/2: Correction in alignment and getting pixel scale
Date: Tue, 17 Jan 2017 23:04:32 +0000 (UTC)

branch: master
commit 4ddc277984260996b4b788fae470372cfca2a17d
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Correction in alignment and getting pixel scale
    
    The conditionals to get the WCS rotation matrix in ImageWarp's
    `makealignmatrix' function were missing an `else' before the `altlin |= 2'
    check! Thus when both flags where on, both conditions were executed which
    would lead the to input matrix begin reset to zero. Therefore no alignment
    would be done.
    
    The reading of WCSLIB's `wcsprm' structure into a proper matrix is an issue
    that is needed in other contexts too. Hence, the new function
    `gal_wcs_array_from_wcsprm' was defined in `lib/wcs.c' to be a central
    place for all Gnuastro to get this matrix when necessary and avoid such
    bugs in the future.
    
    This issue was submitted by Lee Kelvin as bug #50072 along with a FITS
    image as an example. After this correction, the image would be rotated back
    to its original aligned state. However, I noticed that the pixel scales had
    changed! This would not happen for the FITS images that are produced by
    running `make check'!
    
    After going over each step, I noticed that the problem was in the
    `gal_wcs_pixel_scale_deg' function, the pixel scale it returned was wrong
    (even for the aligned image). Until now, the procedure to get the pixel
    scale was very crude: we would use WCSLIB to get the coordinates of the
    points (0,0), (0,1), and (1,0), then use `gal_wcs_angular_distance_deg' to
    find the angular distance between the points. Since a very large number of
    floating point calculations are necessary in the process, I guess it was
    because of floating point errors that we got the wrong pixel scale for
    Lee's image, but correct values for the others.
    
    So, after reviewing my Linear Algebra and also some searches on the
    internet, I found the robust solution to this problem: using `Singular
    Value Decomposition'. Fortunately, the GNU Scientific Library had
    implementations of the decomposition. Thus, as a side-effect to the bug
    reported by Lee, `gal_wcs_pixel_scale_deg' now also returns the correct
    pixel scale through a robust process.
    
    This fixes bug #50072.
---
 THANKS                |    1 +
 bin/imgwarp/imgwarp.c |    9 ++--
 bin/imgwarp/ui.c      |   40 +++++-----------
 lib/gnuastro/wcs.h    |    6 ++-
 lib/wcs.c             |  125 ++++++++++++++++++++++++++++++++++++++++++-------
 5 files changed, 130 insertions(+), 51 deletions(-)

diff --git a/THANKS b/THANKS
index 36d3ad5..036ae65 100644
--- a/THANKS
+++ b/THANKS
@@ -21,6 +21,7 @@ support in Gnuastro. The list is ordered alphabetically.
     Antonio Diaz Diaz                    address@hidden
     Takashi Ichikawa                     address@hidden
     Brandon Invergo                      address@hidden
+    Lee Kelvin                           address@hidden
     Mohammad-Reza Khellat                address@hidden
     Alan Lefor                           address@hidden
     Francesco Montanari                  address@hidden
diff --git a/bin/imgwarp/imgwarp.c b/bin/imgwarp/imgwarp.c
index 1207f58..27c9533 100644
--- a/bin/imgwarp/imgwarp.c
+++ b/bin/imgwarp/imgwarp.c
@@ -430,8 +430,9 @@ correctwcssaveoutput(struct imgwarpparams *p)
 {
   size_t i;
   void *array;
+  double *pixelscale;
+  double *m=p->matrix, diff;
   char keyword[9*FLEN_KEYWORD];
-  double *m=p->matrix, diff, dx, dy;
   struct gal_fits_key_ll *headers=NULL;
   double tpc[4], tcrpix[3], *crpix=p->wcs->crpix, *pc=p->wcs->pc;
   double tinv[4]={p->inverse[0]/p->inverse[8], p->inverse[1]/p->inverse[8],
@@ -484,9 +485,9 @@ correctwcssaveoutput(struct imgwarpparams *p)
      signs are usually different.*/
   if( p->wcs->pc[1]<ABSOLUTEFLTERROR ) p->wcs->pc[1]=0.0f;
   if( p->wcs->pc[2]<ABSOLUTEFLTERROR ) p->wcs->pc[2]=0.0f;
-  gal_wcs_pixel_scale_deg(p->wcs, &dx, &dy);
+  pixelscale=gal_wcs_pixel_scale_deg(p->wcs);
   diff=fabs(p->wcs->pc[0])-fabs(p->wcs->pc[3]);
-  if( fabs(diff/dx)<RELATIVEFLTERROR )
+  if( fabs(diff/pixelscale[0])<RELATIVEFLTERROR )
     p->wcs->pc[3] =  ( (p->wcs->pc[3] < 0.0f ? -1.0f : 1.0f)
                        * fabs(p->wcs->pc[0]) );
 
@@ -495,6 +496,8 @@ correctwcssaveoutput(struct imgwarpparams *p)
                          p->onaxes[1], p->onaxes[0], p->numnul, p->wcs,
                          headers, SPACK_STRING);
 
+  /* Clean up. */
+  free(pixelscale);
   if(array!=p->output)
     free(array);
 }
diff --git a/bin/imgwarp/ui.c b/bin/imgwarp/ui.c
index f8ba672..e3fbebf 100644
--- a/bin/imgwarp/ui.c
+++ b/bin/imgwarp/ui.c
@@ -541,9 +541,7 @@ read_matrix(struct imgwarpparams *p)
 void
 makealignmatrix(struct imgwarpparams *p, double *tmatrix)
 {
-  double A, dx, dy;
-  double amatrix[4];
-  double w[4]={0,0,0,0};
+  double A, *w, *ps, amatrix[4];
 
   /* Check if there is only two WCS axises: */
   if(p->wcs->naxis!=2)
@@ -551,30 +549,12 @@ makealignmatrix(struct imgwarpparams *p, double *tmatrix)
           "axises. For the `--align' option to operate it must be 2",
           p->up.inputname, p->cp.hdu, p->wcs->naxis);
 
-  /* Depending on the type of data, make the input matrix. Note that
-     wcs->altlin is actually bit flags, not integers, so we have to compare
-     with powers of two. */
-  if(p->wcs->altlin |= 1)
-    {
-      w[0]=p->wcs->cdelt[0]*p->wcs->pc[0];
-      w[1]=p->wcs->cdelt[0]*p->wcs->pc[1];
-      w[2]=p->wcs->cdelt[1]*p->wcs->pc[2];
-      w[3]=p->wcs->cdelt[1]*p->wcs->pc[3];
-    }
-  if(p->wcs->altlin |= 2)
-    {
-      w[0]=p->wcs->cd[0];
-      w[1]=p->wcs->cd[1];
-      w[2]=p->wcs->cd[2];
-      w[3]=p->wcs->cd[3];
-    }
-  else
-    error(EXIT_FAILURE, 0, "currently the `--align' option only recognizes "
-          "PCi_ja and CDi_ja keywords, not any others");
 
   /* Find the pixel scale along the two dimensions. Note that we will be
      using the scale along the image X axis for both values. */
-  gal_wcs_pixel_scale_deg(p->wcs, &dx, &dy);
+  w=gal_wcs_array_from_wcsprm(p->wcs);
+  ps=gal_wcs_pixel_scale_deg(p->wcs);
+
 
   /* Lets call the given WCS orientation `W', the rotation matrix we want
      to find as `X' and the final (aligned matrix) to have just one useful
@@ -626,28 +606,32 @@ makealignmatrix(struct imgwarpparams *p, double *tmatrix)
   else
     {
       A = (w[3]/w[1]) - (w[2]/w[0]);
-      amatrix[1] = dx / w[0] / A;
-      amatrix[3] = dx / w[1] / A;
+      amatrix[1] = ps[0] / w[0] / A;
+      amatrix[3] = ps[0] / w[1] / A;
       amatrix[0] = -1 * amatrix[1] * w[3] / w[1];
       amatrix[2] = -1 * amatrix[3] * w[2] / w[0];
     }
 
 
   /* For a check:
-  printf("dx: %e\n", dx);
+  printf("ps: %e\n", ps);
   printf("w:\n");
   printf("  %.8e    %.8e\n", w[0], w[1]);
   printf("  %.8e    %.8e\n", w[2], w[3]);
   printf("x:\n");
   printf("  %.8e    %.8e\n", amatrix[0], amatrix[1]);
   printf("  %.8e    %.8e\n", amatrix[2], amatrix[3]);
+  exit(0);
   */
 
-
   /* Put the matrix elements into the output array: */
   tmatrix[0]=amatrix[0];  tmatrix[1]=amatrix[1]; tmatrix[2]=0.0f;
   tmatrix[3]=amatrix[2];  tmatrix[4]=amatrix[3]; tmatrix[5]=0.0f;
   tmatrix[6]=0.0f;        tmatrix[7]=0.0f;       tmatrix[8]=1.0f;
+
+  /* Clean up. */
+  free(w);
+  free(ps);
 }
 
 
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index aa81d29..4adda93 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -47,6 +47,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 __BEGIN_C_DECLS  /* From C++ preparations */
 
 
+double *
+gal_wcs_array_from_wcsprm(struct wcsprm *wcs);
 
 void
 gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy, double *radec,
@@ -59,8 +61,8 @@ gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double *radec, 
double *xy,
 double
 gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
 
-void
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, double *dy);
+double *
+gal_wcs_pixel_scale_deg(struct wcsprm *wcs);
 
 double
 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
diff --git a/lib/wcs.c b/lib/wcs.c
index 58967bd..2d053d4 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -31,11 +31,66 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <unistd.h>
 #include <assert.h>
 
+#include <gsl/gsl_linalg.h>
+
 #include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 
 
 
+/**************************************************************/
+/**********              Utilities                 ************/
+/**************************************************************/
+double *
+gal_wcs_array_from_wcsprm(struct wcsprm *wcs)
+{
+  double *out;
+  size_t i, j, size=wcs->naxis*wcs->naxis;
+
+  /* Allocate the necessary array. */
+  errno=0;
+  out=malloc(size*sizeof *out);
+  if(out==NULL)
+    error(EXIT_FAILURE, errno, "%zu bytes for `out' in "
+          "`gal_wcs_array_from_wcsprm'", size*sizeof *out);
+
+  /* Fill in the array. */
+  if(wcs->altlin |= 1)          /* Has a PCi_j array. */
+    {
+      for(i=0;i<wcs->naxis;++i)
+        for(j=0;j<wcs->naxis;++j)
+          out[i*wcs->naxis+j] = wcs->cdelt[i] * wcs->pc[i*wcs->naxis+j];
+    }
+  else if(wcs->altlin |= 2)     /* Has CDi_j array */
+    {
+      for(i=0;i<size;++i)
+        out[i]=wcs->cd[i];
+    }
+  else
+    error(EXIT_FAILURE, 0, "currently, `gal_wcs_pixel_scale_deg' only "
+          "recognizes PCi_ja and CDi_ja keywords");
+
+  /* Return the result */
+  return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
 /**************************************************************/
@@ -148,18 +203,42 @@ gal_wcs_angular_distance_deg(double r1, double d1, double 
r2, double d2)
 
 
 /* Return the pixel scale of the image in both dimentions in degrees. */
-void
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, double *dy)
+double *
+gal_wcs_pixel_scale_deg(struct wcsprm *wcs)
 {
-  double radec[6], xy[]={0,0,1,0,0,1};
-
-  /* Get the RA and Dec of the bottom left, bottom right and top left
-     sides of the first pixel in the image. */
-  gal_wcs_xy_array_to_radec(wcs, xy, radec, 3, 2);
-
-  /* Calculate the distances and convert back to degrees: */
-  *dx = gal_wcs_angular_distance_deg(radec[0], radec[1], radec[2], radec[3]);
-  *dy = gal_wcs_angular_distance_deg(radec[0], radec[1], radec[4], radec[5]);
+  gsl_vector S;
+  gsl_matrix A, V;
+  size_t n=wcs->naxis;
+  double *a, *v, *pixscale;
+
+  /* Allocate space for the `v' and `pixscale' arrays. */
+  errno=0; v=malloc(n*n*sizeof *v);
+  if(v==NULL)
+    error(EXIT_FAILURE, errno, "%zu bytes for `v' in "
+          "`gal_wcs_pixel_scale_deg'", n*n*sizeof *v);
+  errno=0; pixscale=malloc(n*sizeof *pixscale);
+  if(pixscale==NULL)
+    error(EXIT_FAILURE, errno, "%zu bytes for `pixscale' in "
+          "`gal_wcs_pixel_scale_deg'", n*sizeof *pixscale);
+
+  /* Write the full matrix into an array, irrespective of what type it was
+     stored in the wcsprm structure (`PCi_j' style or `CDi_j' style). */
+  a=gal_wcs_array_from_wcsprm(wcs);
+
+  /* Fill in the necessary GSL vector and Matrix structures. */
+  S.size=n;     S.stride=1;                S.data=pixscale;
+  V.size1=n;    V.size2=n;    V.tda=n;     V.data=v;
+  A.size1=n;    A.size2=n;    A.tda=n;     A.data=a;
+
+  /* Run GSL's Singular Value Decomposition, using one-sided Jacobi
+     orthogonalization which computes the singular (scale) values to a
+     higher relative accuracy.*/
+  gsl_linalg_SV_decomp_jacobi(&A, &V, &S);
+
+  /* Clean up and return. */
+  free(a);
+  free(v);
+  return pixscale;
 }
 
 
@@ -174,11 +253,21 @@ gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, 
double *dy)
 double
 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
 {
-  double dx, dy;
-
-  /* Get the pixel scales along each axis in degrees. */
-  gal_wcs_pixel_scale_deg(wcs, &dx, &dy);
-
-  /* Return the result */
-  return dx * dy * 3600.0f * 3600.0f;
+  double out;
+  double *pixscale;
+
+  /* A small sanity check. Later, when higher dimensions are necessary, we
+     can find which ones correlate to RA and Dec and use them to find the
+     pixel area in arcsec^2. */
+  if(wcs->naxis!=2)
+    error(EXIT_FAILURE, 0, "`gal_wcs_pixel_area_arcsec2' can currently "
+          "calculate the area only when the image has 2 dimensions.");
+
+  /* Get the pixel scales along each axis in degrees, then multiply. */
+  pixscale=gal_wcs_pixel_scale_deg(wcs);
+
+  /* Clean up and return the result. */
+  out = pixscale[0] * pixscale[1] * 3600.0f * 3600.0f;
+  free(pixscale);
+  return out;
 }



reply via email to

[Prev in Thread] Current Thread [Next in Thread]