[gnuastro-commits] master aae0409 2/2: ConvertType: warning added when images are not aligned

[Mohammad Akhlaghi]

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

    ConvertType: warning added when images are not aligned
    
    Until now, ConvertType would only check the pixel size of the input
    images. However, it could happen that the pixel sizes are the same, but the
    WCS alignment of the images is different.
    
    With this commit, a check has been added at the start of ConvertType for
    this: if the input images have WCS, their central coordinate is checked in
    WCS coordinates and if it different by more than one pixel, ConvertType
    will print an informative warning about it. To provide a fast solution, a
    section has also been added to the ConvertType section of the Gnuastro book
    with a small script that people can just copy-paste to fix the problem.
    
    This issue was found with the help of Raúl Infante Sainz.
---
 bin/convertt/ui.c | 93 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 doc/gnuastro.texi | 66 +++++++++++++++++++++++++++++++++++++--
 2 files changed, 157 insertions(+), 2 deletions(-)

diff --git a/bin/convertt/ui.c b/bin/convertt/ui.c
index 7478aac..219e39e 100644
--- a/bin/convertt/ui.c
+++ b/bin/convertt/ui.c
@@ -627,6 +627,97 @@ ui_make_channels_ll(struct converttparams *p)
 
 
 
+static void
+ui_prepare_input_channels_check_wcs(struct converttparams *p)
+{
+  int printwarning=0;
+  float wcsmatch=1.0;
+  gal_data_t *tmp, *coords=NULL;
+  size_t one=1, numwcs=0, numnonblank=0;
+  double *c1, *c2, r1=NAN, r2=NAN, *pixscale=NULL;
+
+  /* If all the inputs have WCS, check to see if the inputs are aligned and
+     print a warning if they aren't. */
+  for(tmp=p->chll; tmp!=NULL; tmp=tmp->next)
+    {
+      if(tmp->wcs && tmp->type!=GAL_TYPE_INVALID) ++numwcs;
+      if(tmp->type!=GAL_TYPE_INVALID)             ++numnonblank;
+    }
+  if(numwcs==numnonblank)
+    {
+      /* Allocate the coordinate columns. */
+      gal_list_data_add_alloc(&coords, NULL, GAL_TYPE_FLOAT64, 1,
+                              &one, NULL, 0, -1, 1, NULL, NULL, NULL);
+      gal_list_data_add_alloc(&coords, NULL, GAL_TYPE_FLOAT64, 1,
+                              &one, NULL, 0, -1, 1, NULL, NULL, NULL);
+
+      /* Go over each image and put its central pixel in the coordinates
+         and do the world-coordinate transformation. Recall that the C
+         coordinates are the inverse order of FITS coordinates and that
+         FITS coordinates count from 1 (not 0).*/
+      for(tmp=p->chll; tmp!=NULL; tmp=tmp->next)
+        if(tmp->wcs)
+          {
+            /* Fill the coordinate values. */
+            c1=coords->array;
+            c2=coords->next->array;
+            c1[0] = tmp->dsize[1] / 2 + 1;
+            c2[0] = tmp->dsize[0] / 2 + 1;
+
+            /* Get the RA/Dec. */
+            gal_wcs_img_to_world(coords, tmp->wcs, 1);
+
+            /* If the pixel scale hasn't been calculated yet, do it (we
+               only need it once, should be similar in all). */
+            if(pixscale==NULL)
+              pixscale=gal_wcs_pixel_scale(tmp->wcs);
+
+            /* If the reference/first center is not yet defined then write
+               the conversions in it. If it is defined, compare with it
+               with the new dataset and print a warning if necessary. */
+            if( isnan(r1) )
+              { r1=c1[0]; r2=c2[0]; }
+            else
+              {
+                /* For a check.
+                printf("check: %g, %g\n", fabs(c1[0]-r1)/pixscale[0],
+                       fabs(c2[0]-r2)/pixscale[1]);
+                */
+
+                /* See if a warning should be printed. */
+                if( fabs(c1[0]-r1)/pixscale[0] > wcsmatch
+                    || fabs(c2[0]-r2)/pixscale[1] > wcsmatch )
+                  printwarning=1;
+              }
+          }
+    }
+
+  /* Print the warning message if necessary. */
+  if(printwarning && p->cp.quiet==0)
+    {
+      error(EXIT_SUCCESS, 0, "WARNING: The WCS information of the input "
+            "FITS images don't match (by more than %g pixels in the "
+            "center), even though the input images have the same number "
+            "of pixels in each dimension. Therefore the color channels "
+            "of the output colored image may not be aligned. If this is "
+            "not a problem, you can suppress this warning with the "
+            "'--quiet' option.\n\n"
+            "A solution to align your images is provided in the "
+            "\"Aligning images with small WCS offsets\" section of "
+            "Gnuastro's manual. Please run the command below to see "
+            "it (you can return to the command-line by pressing 'q'):\n\n"
+            "   info gnuastro \"Aligning images\"\n",
+            wcsmatch);
+    }
+
+  /* Clean up. */
+  free(pixscale);
+}
+
+
+
+
+
 /* Read the input(s)/channels. */
 static void
 ui_prepare_input_channels(struct converttparams *p)
@@ -702,6 +793,8 @@ ui_prepare_input_channels(struct converttparams *p)
           wcs=tmp->wcs;
       }
 
+  /* Make sure the images are all aligned to the same grid. */
+  ui_prepare_input_channels_check_wcs(p);
 
   /* If ndim is still NULL, then there were no non-blank inputs, so print
      an error. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 3719d6d..323c7c3 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -401,6 +401,7 @@ ConvertType
 
 * Recognized file formats::     Recognized file formats
 * Color::                       Some explanations on color.
+* Aligning images with small WCS offsets::  When the WCS slightly differs.
 * Invoking astconvertt::        Options and arguments to ConvertType.
 
 Table
@@ -9610,6 +9611,7 @@ So before explaining the options and arguments (in 
@ref{Invoking astconvertt}),
 @menu
 * Recognized file formats::     Recognized file formats
 * Color::                       Some explanations on color.
+* Aligning images with small WCS offsets::  When the WCS slightly differs.
 * Invoking astconvertt::        Options and arguments to ConvertType.
 @end menu
 
@@ -9745,7 +9747,7 @@ To print to the standard output, set the output name to 
`@file{stdout}'.
 
 @end table
 
-@node Color, Invoking astconvertt, Recognized file formats, ConvertType
+@node Color, Aligning images with small WCS offsets, Recognized file formats, 
ConvertType
 @subsection Color
 
 @cindex RGB
@@ -9813,7 +9815,67 @@ But thanks to the JPEG compression algorithms, when all 
the pixels of one channe
 Therefore a Grayscale image and a CMYK image that has only the K-channel 
filled are approximately the same file size.
 
 
-@node Invoking astconvertt,  , Color, ConvertType
+@node Aligning images with small WCS offsets, Invoking astconvertt, Color, 
ConvertType
+@subsection Aligning images with small WCS offsets
+
+In order to have nice color images, it is important that the images be 
properly aligned.
+This is usually the case in many scenarios, but it some times happens that the 
images have a small WCS offset, even though they have the same size.
+In such cases you can use the script below to align the images into 
approximately the same pixel grid (to within about 0.5 pixels which is 
sufficient in many color-image usage scenarios).
+
+The script below does the job using Gnuastro's @ref{Warp} and @ref{Crop} 
programs.
+Simply copy the lines below into a plain-text file with your favorite text 
editor and save it as @file{my-align.sh}.
+Don't forget to set the variables of the first three lines to specify the file 
names (without the @file{.fits} suffix) and the HDUs of your inputs.
+These four lines are all you need to edit, leave the rest unchanged.
+Also, if you are copy/pasting the script from a PDF, be careful that the 
single-quotes used in AWK may need to be corrected.
+
+@example
+#!/bin/sh
+
+# Set the input names (without the '.fits' suffix),
+# and their HDUs.
+r=RED_IMAGE_NO_SUFFIX;   rhdu=1
+g=GREEN_IMAGE_NO_SUFFIX; ghdu=1
+b=BLUE_IMAGE_NO_SUFFIX;  bhdu=1
+
+# To stop the script if there is a crash
+set -e
+
+# Align all the images to the celestial poles.
+astwarp $r.fits --align -h$rhdu -o $r-aligned.fits
+astwarp $g.fits --align -h$ghdu -o $g-aligned.fits
+astwarp $b.fits --align -h$bhdu -o $b-aligned.fits
+
+# Calculate the final WCS-based center and image-based width based on
+# the G-band (in RGB) image.
+centerwcs=$(astfits $g-aligned.fits --skycoverage --quiet \
+                    | awk 'NR==1@{printf "%g %g", $1,$2@}')
+widthpix=$(astfits $g-aligned.fits -h1 --quiet \
+                   --keyvalue=NAXIS1,NAXIS2 \
+               | awk '@{printf "%d,%d", $1, $2@}')
+
+# Crop all the images around the desired center and width.
+for f in $r $g $b; do
+  centerpix=$(echo $centerwcs \
+                   | asttable -c'arith $1 $2 wcstoimg' \
+                              --wcsfile=$f-aligned.fits \
+                   | awk '@{printf "%g,%g", $1, $2@}')
+  astcrop $f-aligned.fits --mode=img --width=$widthpix \
+          --center=$centerpix -o$f-use.fits
+  rm $f-aligned.fits
+done
+@end example
+
+Once you have have saved the file and come back to your command-line you can 
run the script like this:
+
+@example
+$ chmod +x my-align.sh
+$ ./my-align.sh
+@end example
+
+@noindent
+Of course, feel free to hack it and modify it to fit your datasets, like the 
rest of Gnuastro, this script is released under GNU GPLv.3 and above, see 
@ref{Your rights}.
+
+@node Invoking astconvertt,  , Aligning images with small WCS offsets, 
ConvertType
 @subsection Invoking ConvertType
 
 ConvertType will convert any recognized input file type to any specified 
output type.