[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 2c554b7 3/4: ImageWarp aligns image and celest
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 2c554b7 3/4: ImageWarp aligns image and celestial axes |
Date: |
Tue, 4 Oct 2016 18:25:15 +0000 (UTC) |
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