[gnuastro-commits] master 2c554b7 3/4: ImageWarp aligns image and celestial axes

[Mohammad Akhlaghi]

branch: master
commit 2c554b778ec581a26fa1d32d6f2db2293ed5a22b
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    ImageWarp aligns image and celestial axes
    
    With the new `--align' (`-a') option added to ImageWarp, it can now align
    the image along the celestial coordinates in the standard fasion in nearly
    all surveys.
    
    This finishes task #14170.
---
 bin/imgwarp/args.h    |   16 ++++----
 bin/imgwarp/imgwarp.c |   15 +++++++-
 bin/imgwarp/imgwarp.h |    9 +++++
 bin/imgwarp/main.h    |    2 +-
 bin/imgwarp/ui.c      |   98 ++++++++++++++++++++++++++++++++++++++++++++++++-
 doc/gnuastro.texi     |   20 +++++++---
 6 files changed, 143 insertions(+), 17 deletions(-)

diff --git a/bin/imgwarp/args.h b/bin/imgwarp/args.h
index 6621e84..e985855 100644
--- a/bin/imgwarp/args.h
+++ b/bin/imgwarp/args.h
@@ -95,6 +95,14 @@ static struct argp_option options[] =
       1
     },
     {
+      "align",
+      'a',
+      0,
+      0,
+      "Align the image and celestial axes.",
+      1
+    },
+    {
       "hstartwcs",
       501,
       "INT",
@@ -110,14 +118,6 @@ static struct argp_option options[] =
       "Header keyword number to stop reading WCS.",
       1
     },
-    {
-      "align",
-      'a',
-      0,
-      0,
-      "Dec on vertical axis and RA on inv. horizontal.",
-      1
-    },
 
 
 
diff --git a/bin/imgwarp/imgwarp.c b/bin/imgwarp/imgwarp.c
index 2786bc3..7be0003 100644
--- a/bin/imgwarp/imgwarp.c
+++ b/bin/imgwarp/imgwarp.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <float.h>
 #include <stdlib.h>
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/polygon.h>
 
@@ -429,8 +430,8 @@ correctwcssaveoutput(struct imgwarpparams *p)
 {
   size_t i;
   void *array;
-  double *m=p->matrix;
   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],
@@ -477,6 +478,18 @@ correctwcssaveoutput(struct imgwarpparams *p)
                                  "element value.", 0, NULL);
     }
 
+  /* Due to floating point errors extremely small values of PC matrix can
+     be set to zero and extremely small differences between PC1_1 and PC2_2
+     can be ignored. The reason for all the `fabs' functions is because the
+     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);
+  diff=fabs(p->wcs->pc[0])-fabs(p->wcs->pc[3]);
+  if( fabs(diff/dx)<RELATIVEFLTERROR )
+    p->wcs->pc[3] =  ( (p->wcs->pc[3] < 0.0f ? -1.0f : 1.0f)
+                       * fabs(p->wcs->pc[0]) );
+
   /* Save the output: */
   gal_fits_array_to_file(p->cp.output, "Warped", p->inputbitpix, array,
                          p->onaxes[1], p->onaxes[0], p->numnul, p->wcs,
diff --git a/bin/imgwarp/imgwarp.h b/bin/imgwarp/imgwarp.h
index c7567b7..5930618 100644
--- a/bin/imgwarp/imgwarp.h
+++ b/bin/imgwarp/imgwarp.h
@@ -25,6 +25,13 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/threads.h>
 
+
+/* Limits to account for floating point errors: */
+#define ABSOLUTEFLTERROR 1e-10
+#define RELATIVEFLTERROR 1e-6
+
+
+/* Internal structure. */
 struct iwpparams
 {
   /* General input parameters: */
@@ -35,6 +42,8 @@ struct iwpparams
   pthread_barrier_t    *b;    /* Barrier to keep threads waiting.      */
 };
 
+
+/* Extenal functions. */
 void
 imgwarp(struct imgwarpparams *p);
 
diff --git a/bin/imgwarp/main.h b/bin/imgwarp/main.h
index c0131f0..3b04619 100644
--- a/bin/imgwarp/main.h
+++ b/bin/imgwarp/main.h
@@ -65,7 +65,7 @@ struct uiparams
 struct imgwarpparams
 {
   /* Other structures: */
-  struct uiparams     up;  /* User interface parameters.                 */
+  struct uiparams      up;  /* User interface parameters.                 */
   struct gal_commonparams cp; /* Common parameters.                      */
 
   /* Input: */
diff --git a/bin/imgwarp/ui.c b/bin/imgwarp/ui.c
index ef4aa24..127d565 100644
--- a/bin/imgwarp/ui.c
+++ b/bin/imgwarp/ui.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <nproc.h>               /* From Gnulib.                   */
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/txtarray.h>
 
@@ -186,6 +187,9 @@ printvalues(FILE *fp, struct imgwarpparams *p)
   if(up->matrixstringset)
     GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("matrix", up->matrixstring);
 
+  if(up->alignset)
+    fprintf(fp, CONF_SHOWFMT"%d\n", "align", up->align);
+
   if(cp->outputset)
     GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("output", cp->output);
 
@@ -304,11 +308,101 @@ readmatrixoption(struct imgwarpparams *p)
 
 
 
-/* Set the matrix so the image is aligned with the axises */
+/* Set the matrix so the image is aligned with the axises. Note that
+   WCSLIB automatically fills the CRPI */
 void
 makealignmatrix(struct imgwarpparams *p)
 {
-  error(EXIT_FAILURE, 0, "The align feature is not implemented yet.");
+  double A, dx, dy;
+  double w[4]={0,0,0,0};
+
+  /* Check if there is only two WCS axises: */
+  if(p->wcs->naxis!=2)
+    error(EXIT_FAILURE, 0, "the WCS structure of %s (hdu: %s) has %d "
+          "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];
+    }
+  else
+    error(EXIT_FAILURE, 0, "currently the `--align' option only recognizes "
+          "PCi_j 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);
+
+  /* Allocate space for the matrix: */
+  errno=0;
+  p->matrix=malloc(4*sizeof *p->matrix);
+  if(p->matrix==NULL)
+    error(EXIT_FAILURE, errno, "%lu bytes for p->matrix",
+          4*sizeof *p->matrix);
+  p->ms0=p->ms1=2;
+
+  /* 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
+     value: `a' (which is the pixel scale):
+
+        x0  x1       w0  w1      -a  0
+        x2  x3   *   w2  w3   =   0  a
+
+     Let's open up the matrix multiplication, so we can find the `X'
+     elements as function of the `W' elements and `a'.
+
+        x0*w0 + x1*w2 = -a                                         (1)
+        x0*w1 + x1*w3 =  0                                         (2)
+        x2*w0 + x3*w2 =  0                                         (3)
+        x2*w1 + x3*w3 =  a                                         (4)
+
+     Let's bring the X with the smaller index in each equation to the left
+     side:
+
+        x0 = (-w2/w0)*x1 - a/w0                                    (5)
+        x0 = (-w3/w1)*x1                                           (6)
+        x2 = (-w2/w0)*x3                                           (7)
+        x2 = (-w3/w1)*x3 + a/w1                                    (8)
+
+    Using (5) and (6) we can find x0 and x1, by first eliminating x0:
+
+       (-w2/w0)*x1 - a/w0 = (-w3/w1)*x1 -> (w3/w1 - w2/w0) * x1 = a/w0
+
+    For easy reading/writing, let's define: A = (w3/w1 - w2/w0)
+
+       --> x1 = a / w0 / A
+       --> x0 = -1 * x1 * w3 / w1
+
+    Similar to the above, we can find x2 and x3 from (7) and (8):
+
+       (-w2/w0)*x3 = (-w3/w1)*x3 + a/w1 -> (w3/w1 - w2/w0) * x3 = a/w1
+
+       --> x3 = a / w1 / A
+       --> x2 = -1 * x3 * w2 / w0
+
+   */
+  A = (w[3]/w[1]) - (w[2]/w[0]);
+  p->matrix[1] = dx / w[0] / A;
+  p->matrix[3] = dx / w[1] / A;
+  p->matrix[0] = -1 * p->matrix[1] * w[3] / w[1];
+  p->matrix[2] = -1 * p->matrix[3] * w[2] / w[0];
+
+  /* For a check:
+  printf("dx: %e\n", dx);
+  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", p->matrix[0], p->matrix[1]);
+  printf("  %.8e    %.8e\n", p->matrix[2], p->matrix[3]);
+  */
 }
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 98fc4a0..c33e2b8 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9139,11 +9139,13 @@ characters. If you want to use the first two, then be 
sure to wrap the
 matrix within double quotation marks (@key{"}) so they are not confused
 with other arguments on the command-line, see @ref{Options}. This also
 applies to values in the configuration files, see @ref{Configuration file
-format}.  The transformation matrix can be either 2 by 2 or 3 by 3
-array. In the former case (if a 2 by 2 matrix is given), then a translation
-is also added to correct for the FITS definition of position to create a
-usable 3 by 3 matrix (see @ref{Warping basics}, @ref{Merging multiple
-warpings}, and the box above for more).
+format}.
+
+The transformation matrix can be either 2 by 2 or 3 by 3 array. In the
+former case (if a 2 by 2 matrix is given), then it is converted to a 3 by 3
+matrix with two translation terms added to correct for the FITS definition
+of position (see @ref{Warping basics}, @ref{Merging multiple warpings}, and
+the box above for more).
 
 @cindex NaN
 The determinant of the matrix has to be non-zero and it must not
@@ -9152,6 +9154,14 @@ elements of the matrix have to be written row by row. So 
for the
 general homography matrix of @ref{Warping basics}, it should be called
 with @command{--matrix=a,b,c,d,e,f,g,h,1}.
 
address@hidden -a
address@hidden --align
+Align the image and celestial (WCS) axes, such that the vertical image
+direction (when viewed in SAO ds9) corresponds to the declination and the
+horizontal axis is the inverse of the Right Ascension (RA). The inverse of
+the RA is chosen so the image can correspond to what you would actually see
+on the sky and is common in most survey images.
+
 @item --hstartwcs
 (@option{=INT}) Specify the first header keyword number (line) that
 should be used to read the WCS information, see the full explanation in